Update 28 Jan 2017

Summary of findings:

  1. % load and total mass are too collinear to include in the same model.
  2. Treatment and % load are too collinear to include in the same model.
  3. For the remainder of the modeling, I’ll use only % load
  4. Conducted all tests for squared terms of continuous predictors (suspect wbf matters)

Differences between model with Treatment vs. model with % load

wingbeat freq = size Met rate = mass carrying amplitude = percent loading

Susie’s (hypotheses): lm( frequency ~ ITspan + trt order)

lm( metabolism ~ total mass)

lm( stroke amplitude ~ %load + trt) (those are probably the same thing)

** Hypotheses: 1. Wingbeat freq moves around for a number of reasons (tired, cold, added mass) 2. Stroke amplitude rescues wingbeat frequency 3. As you get really close to your max, adding more weight doesn’t seem to cost as much.


Analyses of respirometry data:

  1. How does carrying a load during flight affect bumblebees’?
    1. How does load affect respiratory rate?
    2. How does load affect wingbeat frequency?
    3. How does load affect stroke amplitude?
    4. How does load affect wing velocity?
  2. Which flight measurements are most closely associated with respiratory rate in bees – wingbeat frequency, stroke amplitude, or wing velocity?

Setup

Install required packages and read in data Define custom function for evaluating VIF with multilevel models

ipak <- function(pkg){
     new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
     if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
     sapply(pkg, require, character.only = TRUE)
}

packages <- c('ggplot2', 'car', 'lme4', 'gsheet', "MASS", 'influence.ME', 'sjPlot')
ipak(packages)
##      ggplot2          car         lme4       gsheet         MASS 
##         TRUE         TRUE         TRUE         TRUE         TRUE 
## influence.ME       sjPlot 
##         TRUE         TRUE
theme_set(theme_bw())

# read in data -- google sheet called "Bumble mumble grumble"
bdta <- gsheet2tbl("https://docs.google.com/spreadsheets/d/1GUUoFq41Ep3sJNFXiyMcp7fZhYmQCzTcLnhVXdr4WRo/edit?usp=sharing")

Function for calculating vif for lmer

(used later)

vif.mer <- function (fit) {
    ## adapted from rms::vif
    v <- vcov(fit)
    nam <- names(fixef(fit))

    ## exclude intercepts
    ns <- sum(1 * (nam == "Intercept" | nam == "(Intercept)"))
    if (ns > 0) {
        v <- v[-(1:ns), -(1:ns), drop = FALSE]
        nam <- nam[-(1:ns)]
    }
    
    d <- diag(v)^0.5
    v <- diag(solve(v/(d %o% d)))
    names(v) <- nam
    v
}

Data Overview and Engineering

# create centered and squared variables
bdta <- within(bdta, {
     # new variables
     percLoad <- (Mstarved + load) / Mstarved
     totalMass <- Mstarved + load})

bdta <- within(bdta, {
     # centered variables
     percLoad_cent <- scale(bdta$percLoad, center = TRUE, scale = FALSE)[,1]
     totalMass_cent <- scale(totalMass, center = TRUE, scale = FALSE)[,1]
     wbf_cent <- scale(wbf.aud., center = TRUE, scale = FALSE)[,1]
     strokeAmp_cent <- scale(stroke.amplitude, center = TRUE, scale = FALSE)[,1]
     wingVel_cent <- scale(wing.velocity, center = TRUE, scale = FALSE)[,1]
     
     #squared
     percLoad2 <- percLoad_cent^2
     totalMass2 <- totalMass^2
     wbf2 <- wbf_cent^2  
     strokeAmp2 <- strokeAmp_cent^2
     wingVel2 <- wingVel_cent^2
     
     # change factor properties
     # convert trt order to factor
     Treatment.order <- as.factor(as.character(Treatment.order))

     # change reference level to unloaded
     Treatment <- factor(Treatment, levels = c("UL", "L"))
})

summary(bdta)
##     Bee.ID          Treatment.order Treatment    Mstarved     
##  Length:60          1:30            UL:30     Min.   :0.0767  
##  Class :character   2:30            L :30     1st Qu.:0.1213  
##  Mode  :character                             Median :0.1378  
##                                               Mean   :0.1411  
##                                               3rd Qu.:0.1670  
##                                               Max.   :0.1909  
##       load            IT.Span      single.wing.area av.resp..CO2.mL.hr.
##  Min.   :0.01650   Min.   :3.860   Min.   :20.03    Min.   : 5.507     
##  1st Qu.:0.04713   1st Qu.:4.400   1st Qu.:27.93    1st Qu.: 8.935     
##  Median :0.07840   Median :4.620   Median :30.67    Median :10.666     
##  Mean   :0.08332   Mean   :4.594   Mean   :30.70    Mean   :10.600     
##  3rd Qu.:0.11800   3rd Qu.:4.870   3rd Qu.:33.75    3rd Qu.:12.297     
##  Max.   :0.17460   Max.   :5.180   Max.   :39.05    Max.   :16.475     
##     wbf.aud.     stroke.amplitude wing.velocity     totalMass     
##  Min.   :163.9   Min.   : 96.88   Min.   :1.208   Min.   :0.1062  
##  1st Qu.:173.0   1st Qu.:114.16   1st Qu.:1.351   1st Qu.:0.1808  
##  Median :179.1   Median :125.96   Median :1.427   Median :0.2191  
##  Mean   :178.9   Mean   :123.31   Mean   :1.465   Mean   :0.2245  
##  3rd Qu.:185.3   3rd Qu.:133.61   3rd Qu.:1.588   3rd Qu.:0.2646  
##  Max.   :194.7   Max.   :144.78   Max.   :1.843   Max.   :0.3550  
##     percLoad        wingVel2           strokeAmp2      
##  Min.   :1.143   Min.   :7.978e-05   Min.   :  0.0093  
##  1st Qu.:1.321   1st Qu.:3.980e-03   1st Qu.: 42.9577  
##  Median :1.576   Median :1.444e-02   Median : 99.2962  
##  Mean   :1.600   Mean   :2.509e-02   Mean   :148.2959  
##  3rd Qu.:1.863   3rd Qu.:3.750e-02   3rd Qu.:188.1689  
##  Max.   :2.180   Max.   :1.423e-01   Max.   :698.3210  
##       wbf2             totalMass2        percLoad2        
##  Min.   :  0.00127   Min.   :0.01128   Min.   :0.0003255  
##  1st Qu.:  6.59393   1st Qu.:0.03270   1st Qu.:0.0353246  
##  Median : 38.46046   Median :0.04801   Median :0.0764175  
##  Mean   : 57.11835   Mean   :0.05374   Mean   :0.0964447  
##  3rd Qu.: 77.70129   3rd Qu.:0.07003   3rd Qu.:0.1477398  
##  Max.   :250.38161   Max.   :0.12602   Max.   :0.3363085  
##   wingVel_cent      strokeAmp_cent       wbf_cent      
##  Min.   :-0.25706   Min.   :-26.426   Min.   :-15.074  
##  1st Qu.:-0.11484   1st Qu.: -9.145   1st Qu.: -5.972  
##  Median :-0.03848   Median :  2.656   Median :  0.185  
##  Mean   : 0.00000   Mean   :  0.000   Mean   :  0.000  
##  3rd Qu.: 0.12225   3rd Qu.: 10.307   3rd Qu.:  6.394  
##  Max.   : 0.37723   Max.   : 21.473   Max.   : 15.823  
##  totalMass_cent      percLoad_cent     
##  Min.   :-0.118253   Min.   :-0.45679  
##  1st Qu.:-0.043628   1st Qu.:-0.27912  
##  Median :-0.005353   Median :-0.02372  
##  Mean   : 0.000000   Mean   : 0.00000  
##  3rd Qu.: 0.040172   3rd Qu.: 0.26297  
##  Max.   : 0.130547   Max.   : 0.57992
ggplot(bdta, aes(x = percLoad, fill = Treatment)) + 
     geom_histogram()  
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# seems like percent load and trt are very related
# visualize % load vs total mass
ggplot(bdta, aes(x = percLoad, y = totalMass, color = Treatment)) + 
     geom_point()

# seem to be strongly associated
# see number of observations per bee
table(bdta$Bee.ID) # each bee has two observations
## 
## E01 E03 E05 E09 E10 E11 E12 E16 E20 E24 E26 E28 E30 E31 E32 E33 E35 E36 
##   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2 
## E37 E39 E41 E42 E44 E45 E49 E50 E53 E54 E56 E59 
##   2   2   2   2   2   2   2   2   2   2   2   2
# get number of unque bees
length(table(bdta$Bee.ID))
## [1] 30
# get treatment orders
trtOrders_loadedSecond <- sapply(X = unique(bdta$Bee.ID), FUN = function(x){
     tmp = bdta[bdta$Bee.ID == x, c("Treatment.order", "Treatment")]
     if("2_L" %in% paste(tmp[,1], tmp[,2], sep = "_")){
          loadedSecond = TRUE
     }
     else loadedSecond = FALSE
     return(loadedSecond)
})

table(trtOrders_loadedSecond)
## trtOrders_loadedSecond
## FALSE 
##    30
# visualize scatterplot
car::scatterplotMatrix(bdta[, 4:13])

Variance Inflation Factor and PCA

# VIF  is high among the bee size predictors
# note: this function is almost the same as car::vif
vif.mer(lmer(av.resp..CO2.mL.hr. ~ Treatment + Treatment.order + Mstarved + IT.Span + totalMass +  single.wing.area + percLoad +  (1|Bee.ID), data = bdta))
##       TreatmentL Treatment.order2         Mstarved          IT.Span 
##        16.591645         1.270430        23.112962         6.658432 
##        totalMass single.wing.area         percLoad 
##        32.931859        10.903672        38.415618
car::scatterplotMatrix(bdta[, c("Mstarved", "IT.Span", "single.wing.area")])

# principle components
aa = prcomp(bdta[, c("Mstarved", "IT.Span", "single.wing.area")], center = TRUE, scale = TRUE)
summary(aa) # 1st pc explains ~95% of the variance in the three measurements of size
## Importance of components:
##                           PC1     PC2     PC3
## Standard deviation     1.6852 0.33843 0.21323
## Proportion of Variance 0.9467 0.03818 0.01516
## Cumulative Proportion  0.9467 0.98484 1.00000
biplot(aa) # shows that all three size measurement are correlated

# note, I changed the signs of the predictions so that higher PC1 values 
# correspond to bigger bees
p1 = -predict(aa)[,1] 

# add PC1 scores to dataset
bdta$size_pc1 = p1

# show scatterplot matrix to see correlations among size predictors
car::scatterplotMatrix(bdta[, c("Mstarved", "IT.Span", "single.wing.area", "size_pc1", "percLoad")])

# check VIF one more time
# VIF  is high among the bee size predictors
vif.mer(lmer(av.resp..CO2.mL.hr. ~ Treatment + Treatment.order + size_pc1 + percLoad + totalMass +  (1|Bee.ID), data = bdta))
##       TreatmentL Treatment.order2         size_pc1         percLoad 
##        16.643643         1.258990         7.448179        30.631628 
##        totalMass 
##        26.790114
# remove treatment 
vif.mer(lmer(av.resp..CO2.mL.hr. ~ Treatment.order + size_pc1 + percLoad + totalMass + (1|Bee.ID), data = bdta))
## Treatment.order2         size_pc1         percLoad        totalMass 
##         1.014834         7.698050        20.910198        25.813622
# looks like we cannot have total mass and percentLoad in the same model, according to VIF

# remove total mass 
vif.mer(lmer(av.resp..CO2.mL.hr. ~ Treatment.order + size_pc1 + percLoad +  (1|Bee.ID), data = bdta))
## Treatment.order2         size_pc1         percLoad 
##         1.003811         1.004814         1.008624

Multilevel models


Resp Rate

# make a full model with all two-way interactions
m1 <- lmer(av.resp..CO2.mL.hr. ~  (size_pc1 +  Treatment.order + percLoad_cent)^2 + 
                percLoad2 + 
           + (1|Bee.ID), data = bdta)

summary(m1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## av.resp..CO2.mL.hr. ~ (size_pc1 + Treatment.order + percLoad_cent)^2 +  
##     percLoad2 + +(1 | Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: 182.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7335 -0.4895  0.1150  0.4689  2.3779 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 1.4520   1.2050  
##  Residual             0.5021   0.7086  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                Estimate Std. Error t value
## (Intercept)                    11.23228    0.27961   40.17
## size_pc1                        1.28130    0.15748    8.14
## Treatment.order2               -0.28829    0.18935   -1.52
## percLoad_cent                   4.90035    0.70897    6.91
## percLoad2                      -5.10731    1.38346   -3.69
## size_pc1:Treatment.order2      -0.17569    0.12108   -1.45
## size_pc1:percLoad_cent          0.09533    0.21461    0.44
## Treatment.order2:percLoad_cent -1.19963    1.15631   -1.04
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2 s_1:T. s_1:L_
## size_pc1    -0.001                                          
## Trtmnt.rdr2 -0.253  0.032                                   
## percLod_cnt -0.042  0.233  0.066                            
## percLoad2   -0.392 -0.049 -0.204 -0.043                     
## sz_pc1:Tr.2 -0.049 -0.410 -0.089 -0.325  0.219              
## sz_pc1:prL_ -0.030 -0.105 -0.205 -0.146  0.331  0.289       
## Trtmnt.2:L_  0.071 -0.207 -0.029 -0.886 -0.035  0.284  0.159
### LRT's for interactions
m2.0 <- update(m1, .~. - size_pc1:percLoad_cent)

anova(m1, m2.0) # likelihood ratio test for interaction of treatment order
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.0: av.resp..CO2.mL.hr. ~ size_pc1 + Treatment.order + percLoad_cent + 
## m2.0:     percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order + Treatment.order:percLoad_cent
## m1: av.resp..CO2.mL.hr. ~ (size_pc1 + Treatment.order + percLoad_cent)^2 + 
## m1:     percLoad2 + +(1 | Bee.ID)
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m2.0  9 195.92 214.77 -88.961   177.92                         
## m1   10 197.68 218.62 -88.841   177.68 0.2397      1     0.6244
summary(m2.0)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## av.resp..CO2.mL.hr. ~ size_pc1 + Treatment.order + percLoad_cent +  
##     percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order + Treatment.order:percLoad_cent
##    Data: bdta
## 
## REML criterion at convergence: 181.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.72009 -0.52773  0.07283  0.48826  2.35677 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 1.4603   1.2084  
##  Residual             0.4862   0.6973  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                Estimate Std. Error t value
## (Intercept)                     11.2371     0.2783   40.37
## size_pc1                         1.2881     0.1562    8.24
## Treatment.order2                -0.2709     0.1824   -1.49
## percLoad_cent                    4.9339     0.6939    7.11
## percLoad2                       -5.3211     1.2854   -4.14
## size_pc1:Treatment.order2       -0.1907     0.1141   -1.67
## Treatment.order2:percLoad_cent  -1.2596     1.1316   -1.11
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2 s_1:T.
## size_pc1    -0.004                                   
## Trtmnt.rdr2 -0.262  0.010                            
## percLod_cnt -0.048  0.219  0.036                     
## percLoad2   -0.401 -0.014 -0.147  0.008              
## sz_pc1:Tr.2 -0.041 -0.393 -0.031 -0.299  0.136       
## Trtmnt.2:L_  0.078 -0.192  0.004 -0.885 -0.096  0.253
m2.1 <- update(m2.0,.~. - Treatment.order:percLoad_cent)
anova(m2.0, m2.1)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.1: av.resp..CO2.mL.hr. ~ size_pc1 + Treatment.order + percLoad_cent + 
## m2.1:     percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order
## m2.0: av.resp..CO2.mL.hr. ~ size_pc1 + Treatment.order + percLoad_cent + 
## m2.0:     percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order + Treatment.order:percLoad_cent
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m2.1  8 195.15 211.90 -89.573   179.15                         
## m2.0  9 195.92 214.77 -88.961   177.92 1.2255      1     0.2683
summary(m2.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## av.resp..CO2.mL.hr. ~ size_pc1 + Treatment.order + percLoad_cent +  
##     percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order
##    Data: bdta
## 
## REML criterion at convergence: 185
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.69131 -0.46047  0.00547  0.49254  2.34469 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 1.5578   1.2481  
##  Residual             0.4642   0.6813  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                           Estimate Std. Error t value
## (Intercept)                11.2638     0.2811   40.07
## size_pc1                    1.2546     0.1562    8.03
## Treatment.order2           -0.2696     0.1782   -1.51
## percLoad_cent               4.2487     0.3162   13.44
## percLoad2                  -5.4887     1.2518   -4.38
## size_pc1:Treatment.order2  -0.1589     0.1079   -1.47
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2
## size_pc1     0.010                            
## Trtmnt.rdr2 -0.254  0.011                     
## percLod_cnt  0.044  0.102  0.086              
## percLoad2   -0.383 -0.032 -0.147 -0.165       
## sz_pc1:Tr.2 -0.061 -0.348 -0.033 -0.167  0.167
m2.2 <- update(m2.1, .~. - size_pc1:Treatment.order)
anova(m2.1, m2.2) ## drop all interactions
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.2: av.resp..CO2.mL.hr. ~ size_pc1 + Treatment.order + percLoad_cent + 
## m2.2:     percLoad2 + (1 | Bee.ID)
## m2.1: av.resp..CO2.mL.hr. ~ size_pc1 + Treatment.order + percLoad_cent + 
## m2.1:     percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m2.2  7 195.54 210.21 -90.773   181.54                         
## m2.1  8 195.15 211.90 -89.573   179.15 2.3985      1     0.1215
summary(m2.2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## av.resp..CO2.mL.hr. ~ size_pc1 + Treatment.order + percLoad_cent +  
##     percLoad2 + (1 | Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: 184.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6892 -0.4659  0.0722  0.3830  2.1677 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 1.5420   1.2418  
##  Residual             0.4856   0.6968  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)       11.2373     0.2818   39.88
## size_pc1           1.1746     0.1462    8.03
## Treatment.order2  -0.2786     0.1822   -1.53
## percLoad_cent      4.1717     0.3187   13.09
## percLoad2         -5.1672     1.2616   -4.10
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2 prcLd_
## size_pc1    -0.012                     
## Trtmnt.rdr2 -0.261 -0.001              
## percLod_cnt  0.035  0.049  0.081       
## percLoad2   -0.385  0.029 -0.144 -0.141
# renaming model to simplify later typing
m2 <- m2.2

##### LRTs for main effects
## Treatment Order
m3 <- update(m2, .~. - Treatment.order)
anova(m2, m3, test = "Chi") # drop trt order (different than model with treatment)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m3: av.resp..CO2.mL.hr. ~ size_pc1 + percLoad_cent + percLoad2 + 
## m3:     (1 | Bee.ID)
## m2: av.resp..CO2.mL.hr. ~ size_pc1 + Treatment.order + percLoad_cent + 
## m2:     percLoad2 + (1 | Bee.ID)
##    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m3  6 196.04 208.60 -92.019   184.04                         
## m2  7 195.54 210.21 -90.773   181.54 2.4925      1     0.1144
summary(m3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ size_pc1 + percLoad_cent + percLoad2 +  
##     (1 | Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: 185.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.61825 -0.42513 -0.04217  0.52054  2.29988 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 1.5364   1.2395  
##  Residual             0.5073   0.7123  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept)    11.1235     0.2735   40.67
## size_pc1        1.1745     0.1465    8.02
## percLoad_cent   4.2122     0.3245   12.98
## percLoad2      -5.4323     1.2754   -4.26
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 prcLd_
## size_pc1    -0.013              
## percLod_cnt  0.059  0.050       
## percLoad2   -0.450  0.029 -0.131
# LRT for size
m4 <- update(m3, .~. - percLoad2)
anova(m3, m4) # keep percLoad2
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m4: av.resp..CO2.mL.hr. ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
## m3: av.resp..CO2.mL.hr. ~ size_pc1 + percLoad_cent + percLoad2 + 
## m3:     (1 | Bee.ID)
##    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## m4  5 208.25 218.73 -99.127   198.25                             
## m3  6 196.04 208.60 -92.019   184.04 14.216      1  0.0001629 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ size_pc1 + percLoad_cent + percLoad2 +  
##     (1 | Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: 185.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.61825 -0.42513 -0.04217  0.52054  2.29988 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 1.5364   1.2395  
##  Residual             0.5073   0.7123  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept)    11.1235     0.2735   40.67
## size_pc1        1.1745     0.1465    8.02
## percLoad_cent   4.2122     0.3245   12.98
## percLoad2      -5.4323     1.2754   -4.26
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 prcLd_
## size_pc1    -0.013              
## percLod_cnt  0.059  0.050       
## percLoad2   -0.450  0.029 -0.131
# LRT for Treatment (load)
m5 <- update(m3, .~. - size_pc1)
anova(m3, m5) # keep size
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m5: av.resp..CO2.mL.hr. ~ percLoad_cent + percLoad2 + (1 | Bee.ID)
## m3: av.resp..CO2.mL.hr. ~ size_pc1 + percLoad_cent + percLoad2 + 
## m3:     (1 | Bee.ID)
##    Df    AIC    BIC   logLik deviance  Chisq Chi Df Pr(>Chisq)    
## m5  5 229.84 240.31 -109.920   219.84                             
## m3  6 196.04 208.60  -92.019   184.04 35.801      1  2.185e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m6 <- update(m3, .~. - percLoad_cent)
anova(m3, m6) # keep squared term
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m6: av.resp..CO2.mL.hr. ~ size_pc1 + percLoad2 + (1 | Bee.ID)
## m3: av.resp..CO2.mL.hr. ~ size_pc1 + percLoad_cent + percLoad2 + 
## m3:     (1 | Bee.ID)
##    Df    AIC    BIC   logLik deviance  Chisq Chi Df Pr(>Chisq)    
## m6  5 257.02 267.49 -123.510   247.02                             
## m3  6 196.04 208.60  -92.019   184.04 62.982      1  2.086e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# summarize final model for paper
summary(m3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ size_pc1 + percLoad_cent + percLoad2 +  
##     (1 | Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: 185.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.61825 -0.42513 -0.04217  0.52054  2.29988 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 1.5364   1.2395  
##  Residual             0.5073   0.7123  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept)    11.1235     0.2735   40.67
## size_pc1        1.1745     0.1465    8.02
## percLoad_cent   4.2122     0.3245   12.98
## percLoad2      -5.4323     1.2754   -4.26
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 prcLd_
## size_pc1    -0.013              
## percLod_cnt  0.059  0.050       
## percLoad2   -0.450  0.029 -0.131
# write output
summary(m3)$coefficients  
##                Estimate Std. Error   t value
## (Intercept)   11.123542  0.2734947 40.671875
## size_pc1       1.174509  0.1464512  8.019799
## percLoad_cent  4.212203  0.3245350 12.979195
## percLoad2     -5.432315  1.2754342 -4.259189
write.csv( summary(m3)$coefficients, file = "RespCoefs_percLoad.csv" )

resp rate diagnostics

# qq plot
qqnorm(resid(m3), main = "")
qqline(resid(m3)) # good

# residual plot
plot(fitted(m3), resid(m3), xlab = "fitted", ylab = "residuals")
abline(0,0)

# check cook's distance -- should be less than 1, so we're good
infl <- influence(m3, obs = TRUE)
plot(infl, which = 'cook')

# QQPlot for group-level effects
qqnorm(ranef(m3)$Bee.ID[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(m3)$Bee.ID[[1]]) # looks good

resp rate visualization

newdf <- data.frame(size_pc1 = 0, Treatment.order = factor(1), percLoad_cent = seq(min(bdta$percLoad_cent), max(bdta$percLoad_cent), length = 100))
newdf$percLoad2 <- newdf$percLoad_cent^2

preds1 <- predict(m3, newdf, re.form = NA)

# plot of prediction for average sized bee (line)
# raw data plotted as points
ggplot(bdta, aes(x = percLoad_cent + mean(bdta$percLoad) , y = av.resp..CO2.mL.hr.)) + 
     geom_point(alpha = 0.3) + 
     geom_line(data = newdf, aes(x = percLoad_cent + mean(bdta$percLoad), y = preds1)) + 
     labs(x = "Mass during flight (% of starved mass in g)")

Wingbeat Freq

# make a full model with all two-way interactions
m1 <- lmer(wbf.aud. ~  (size_pc1 +  Treatment.order + percLoad_cent)^2 + 
                percLoad2  + 
                (1|Bee.ID), data = bdta)

summary(m1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## wbf.aud. ~ (size_pc1 + Treatment.order + percLoad_cent)^2 + percLoad2 +  
##     (1 | Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: 326.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8304 -0.7064  0.1416  0.5616  1.5135 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 17.295   4.159   
##  Residual              9.876   3.143   
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                Estimate Std. Error t value
## (Intercept)                    183.2406     1.0745  170.53
## size_pc1                        -2.9657     0.5916   -5.01
## Treatment.order2                -2.6627     0.8394   -3.17
## percLoad_cent                    7.0660     2.9186    2.42
## percLoad2                      -31.4157     6.0747   -5.17
## size_pc1:Treatment.order2        0.2089     0.5339    0.39
## size_pc1:percLoad_cent          -0.5715     0.9461   -0.60
## Treatment.order2:percLoad_cent  -0.7048     4.6314   -0.15
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2 s_1:T. s_1:L_
## size_pc1     0.005                                          
## Trtmnt.rdr2 -0.294  0.038                                   
## percLod_cnt -0.026  0.256  0.072                            
## percLoad2   -0.446 -0.066 -0.202 -0.081                     
## sz_pc1:Tr.2 -0.062 -0.476 -0.089 -0.311  0.228              
## sz_pc1:prL_ -0.035 -0.121 -0.204 -0.146  0.332  0.286       
## Trtmnt.2:L_  0.061 -0.224 -0.034 -0.868  0.001  0.267  0.161
### LRT's for interactions
m2.0 <- update(m1, .~. - Treatment.order:percLoad_cent)

anova(m1, m2.0)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.0: wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 + 
## m2.0:     (1 | Bee.ID) + size_pc1:Treatment.order + size_pc1:percLoad_cent
## m1: wbf.aud. ~ (size_pc1 + Treatment.order + percLoad_cent)^2 + percLoad2 + 
## m1:     (1 | Bee.ID)
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m2.0  9 362.18 381.03 -172.09   344.18                         
## m1   10 364.15 385.10 -172.08   344.15 0.0311      1       0.86
summary(m2.0)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 +  
##     (1 | Bee.ID) + size_pc1:Treatment.order + size_pc1:percLoad_cent
##    Data: bdta
## 
## REML criterion at convergence: 331.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8513 -0.6994  0.1341  0.5629  1.5305 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 16.921   4.113   
##  Residual              9.716   3.117   
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                           Estimate Std. Error t value
## (Intercept)               183.2499     1.0623  172.50
## size_pc1                   -2.9860     0.5709   -5.23
## Treatment.order2           -2.6673     0.8321   -3.21
## percLoad_cent               6.6792     1.4349    4.65
## percLoad2                 -31.4059     6.0246   -5.21
## size_pc1:Treatment.order2   0.2308     0.5103    0.45
## size_pc1:percLoad_cent     -0.5481     0.9261   -0.59
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2 s_1:T.
## size_pc1     0.020                                   
## Trtmnt.rdr2 -0.293  0.031                            
## percLod_cnt  0.054  0.128  0.086                     
## percLoad2   -0.447 -0.068 -0.202 -0.162              
## sz_pc1:Tr.2 -0.082 -0.444 -0.082 -0.164  0.237       
## sz_pc1:prL_ -0.045 -0.088 -0.201 -0.012  0.336  0.255
m2.1 <- update(m2.0,.~. - size_pc1:Treatment.order)
anova(m2.0, m2.1)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.1: wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 + 
## m2.1:     (1 | Bee.ID) + size_pc1:percLoad_cent
## m2.0: wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 + 
## m2.0:     (1 | Bee.ID) + size_pc1:Treatment.order + size_pc1:percLoad_cent
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m2.1  8 360.42 377.17 -172.21   344.42                         
## m2.0  9 362.18 381.03 -172.09   344.18 0.2349      1     0.6279
summary(m2.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 +  
##     (1 | Bee.ID) + size_pc1:percLoad_cent
##    Data: bdta
## 
## REML criterion at convergence: 332.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8110 -0.6782  0.1540  0.5820  1.4972 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 17.141   4.140   
##  Residual              9.392   3.065   
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                        Estimate Std. Error t value
## (Intercept)            183.2949     1.0536  173.96
## size_pc1                -2.8715     0.5122   -5.61
## Treatment.order2        -2.6345     0.8154   -3.23
## percLoad_cent            6.7941     1.3926    4.88
## percLoad2              -32.1204     5.7604   -5.58
## size_pc1:percLoad_cent  -0.6559     0.8809   -0.74
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2
## size_pc1    -0.018                            
## Trtmnt.rdr2 -0.298 -0.006                     
## percLod_cnt  0.041  0.062  0.074              
## percLoad2   -0.437  0.042 -0.189 -0.128       
## sz_pc1:prL_ -0.025  0.028 -0.187  0.031  0.294
m2.2 <- update(m2.1, .~. - size_pc1:percLoad_cent)
anova(m2.1, m2.2) ## drop all interactions
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.2: wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 + 
## m2.2:     (1 | Bee.ID)
## m2.1: wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 + 
## m2.1:     (1 | Bee.ID) + size_pc1:percLoad_cent
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m2.2  7 359.05 373.71 -172.53   345.05                         
## m2.1  8 360.42 377.17 -172.21   344.42 0.6319      1     0.4267
summary(m2.2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 +  
##     (1 | Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: 334.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8541 -0.7036  0.1281  0.5915  1.3890 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 17.163   4.143   
##  Residual              9.258   3.043   
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)      183.2771     1.0500  174.54
## size_pc1          -2.8608     0.5114   -5.59
## Treatment.order2  -2.7473     0.7953   -3.45
## percLoad_cent      6.8290     1.3823    4.94
## percLoad2        -30.8797     5.4673   -5.65
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2 prcLd_
## size_pc1    -0.017                     
## Trtmnt.rdr2 -0.307 -0.001              
## percLod_cnt  0.041  0.061  0.081       
## percLoad2   -0.448  0.035 -0.143 -0.144
# renaming model to simplify later typing
m2 <- m2.2

##### LRTs for main effects
## Treatment Order
m3 <- update(m2, .~. - Treatment.order)
anova(m2, m3, test = "Chi") # keep treatment order
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m3: wbf.aud. ~ size_pc1 + percLoad_cent + percLoad2 + (1 | Bee.ID)
## m2: wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 + 
## m2:     (1 | Bee.ID)
##    Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)    
## m3  6 368.10 380.67 -178.05   356.10                            
## m2  7 359.05 373.71 -172.53   345.05 11.05      1  0.0008867 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# LRT for size
m4 <- update(m2, .~. - size_pc1)
anova(m2, m4) # keep size_pc1
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m4: wbf.aud. ~ Treatment.order + percLoad_cent + percLoad2 + (1 | 
## m4:     Bee.ID)
## m2: wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 + 
## m2:     (1 | Bee.ID)
##    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## m4  6 379.62 392.19 -183.81   367.62                             
## m2  7 359.05 373.71 -172.53   345.05 22.571      1  2.025e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# LRT for Treatment (load)
m5 <- update(m2, .~. - percLoad_cent)
anova(m2, m5) # keep percLoad
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m5: wbf.aud. ~ size_pc1 + Treatment.order + percLoad2 + (1 | Bee.ID)
## m2: wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 + 
## m2:     (1 | Bee.ID)
##    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## m5  6 375.91 388.47 -181.95   363.91                             
## m2  7 359.05 373.71 -172.53   345.05 18.858      1  1.408e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m6 <- update(m2, .~. - percLoad2)
anova(m2, m6) # keep squared term
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m6: wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + (1 | 
## m6:     Bee.ID)
## m2: wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 + 
## m2:     (1 | Bee.ID)
##    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## m6  6 378.98 391.54 -183.49   366.98                             
## m2  7 359.05 373.71 -172.53   345.05 21.927      1  2.832e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# summarize final model for paper
summary(m2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 +  
##     (1 | Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: 334.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8541 -0.7036  0.1281  0.5915  1.3890 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 17.163   4.143   
##  Residual              9.258   3.043   
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)      183.2771     1.0500  174.54
## size_pc1          -2.8608     0.5114   -5.59
## Treatment.order2  -2.7473     0.7953   -3.45
## percLoad_cent      6.8290     1.3823    4.94
## percLoad2        -30.8797     5.4673   -5.65
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2 prcLd_
## size_pc1    -0.017                     
## Trtmnt.rdr2 -0.307 -0.001              
## percLod_cnt  0.041  0.061  0.081       
## percLoad2   -0.448  0.035 -0.143 -0.144
# write output
summary(m2)$coefficients  
##                    Estimate Std. Error    t value
## (Intercept)      183.277063  1.0500400 174.542933
## size_pc1          -2.860783  0.5114480  -5.593497
## Treatment.order2  -2.747317  0.7952937  -3.454469
## percLoad_cent      6.829045  1.3823068   4.940325
## percLoad2        -30.879690  5.4672717  -5.648099
write.csv( summary(m2)$coefficients, file = "FreqCoefs_percLoad.csv" )

wbf diagnostics

# qq plot
qqnorm(resid(m2), main = "")
qqline(resid(m2)) # good

# check influence
plot(influence(m2, obs = TRUE ), which = 'cook')

# residual plot
plot(fitted(m2), residuals(m2, "deviance"), xlab = "fitted", ylab = "residuals")
abline(0,0)

# QQPlot for group-level effects
qqnorm(ranef(m2)$Bee.ID[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(m2)$Bee.ID[[1]]) # looks good

wbf visualization

newdf <- data.frame(size_pc1 = 0, Treatment.order = factor(1), percLoad_cent = seq(min(bdta$percLoad_cent), max(bdta$percLoad_cent), length = 100))
newdf$percLoad2 <- newdf$percLoad_cent^2

preds1 <- predict(m2, newdf, re.form = NA) #predict for mean of random effect

# plot of prediction for average sized bee (line)
# raw data plotted as points
ggplot(bdta, aes(x = percLoad_cent + mean(bdta$percLoad) , y = wbf.aud.)) + 
     geom_point(alpha = 0.3) + 
     geom_line(data = newdf, aes(x = percLoad_cent + mean(bdta$percLoad), y = preds1)) + 
     labs(x = "Mass during flight (% of starved mass in g)")

Stroke Amplitude

# make a full model with all two-way interactions
m1 <- lmer(stroke.amplitude ~  (size_pc1 +  Treatment.order + percLoad_cent)^2 + 
                percLoad2 + 
                (1|Bee.ID), data = bdta)

summary(m1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## stroke.amplitude ~ (size_pc1 + Treatment.order + percLoad_cent)^2 +  
##     percLoad2 + (1 | Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: 352.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.09359 -0.52135  0.01418  0.41271  2.08441 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 13.01    3.606   
##  Residual             24.48    4.948   
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                Estimate Std. Error t value
## (Intercept)                    124.1538     1.3504   91.94
## size_pc1                         1.9451     0.7042    2.76
## Treatment.order2                -0.4425     1.3197   -0.34
## percLoad_cent                   39.0584     3.8975   10.02
## percLoad2                       -9.2460     9.3003   -0.99
## size_pc1:Treatment.order2       -1.3443     0.8311   -1.62
## size_pc1:percLoad_cent          -3.1048     1.4603   -2.13
## Treatment.order2:percLoad_cent  -5.6584     5.7085   -0.99
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2 s_1:T. s_1:L_
## size_pc1     0.025                                          
## Trtmnt.rdr2 -0.374  0.049                                   
## percLod_cnt  0.026  0.292  0.087                            
## percLoad2   -0.540 -0.105 -0.198 -0.167                     
## sz_pc1:Tr.2 -0.088 -0.611 -0.087 -0.284  0.242              
## sz_pc1:prL_ -0.040 -0.151 -0.199 -0.147  0.329  0.278       
## Trtmnt.2:L_  0.016 -0.244 -0.047 -0.821  0.086  0.237  0.172
### LRT's for interactions
m2.0 <- update(m1, .~. - Treatment.order:percLoad_cent )

anova(m1, m2.0)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.0: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent + 
## m2.0:     percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order + size_pc1:percLoad_cent
## m1: stroke.amplitude ~ (size_pc1 + Treatment.order + percLoad_cent)^2 + 
## m1:     percLoad2 + (1 | Bee.ID)
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m2.0  9 394.21 413.05 -188.10   376.21                         
## m1   10 395.22 416.16 -187.61   375.22 0.9852      1     0.3209
summary(m2.0)
## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent +  
##     percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order + size_pc1:percLoad_cent
##    Data: bdta
## 
## REML criterion at convergence: 358.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.06103 -0.56098 -0.01219  0.50169  2.15159 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 14.18    3.765   
##  Residual             23.65    4.864   
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                           Estimate Std. Error t value
## (Intercept)               124.1660     1.3486   92.07
## size_pc1                    1.7778     0.6856    2.59
## Treatment.order2           -0.5058     1.2961   -0.39
## percLoad_cent              35.9655     2.1941   16.39
## percLoad2                  -8.3350     9.1404   -0.91
## size_pc1:Treatment.order2  -1.1503     0.7940   -1.45
## size_pc1:percLoad_cent     -2.8411     1.4174   -2.00
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2 s_1:T.
## size_pc1     0.029                                   
## Trtmnt.rdr2 -0.367  0.038                            
## percLod_cnt  0.068  0.163  0.085                     
## percLoad2   -0.536 -0.086 -0.196 -0.168              
## sz_pc1:Tr.2 -0.094 -0.575 -0.078 -0.162  0.229       
## sz_pc1:prL_ -0.044 -0.112 -0.195 -0.011  0.322  0.248
m2.1 <- update(m2.0,.~. - size_pc1:Treatment.order)
anova(m2.0, m2.1)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.1: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent + 
## m2.1:     percLoad2 + (1 | Bee.ID) + size_pc1:percLoad_cent
## m2.0: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent + 
## m2.0:     percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order + size_pc1:percLoad_cent
##      Df    AIC    BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## m2.1  8 394.61 411.36 -189.3   378.61                         
## m2.0  9 394.21 413.05 -188.1   376.21 2.4002      1     0.1213
summary(m2.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent +  
##     percLoad2 + (1 | Bee.ID) + size_pc1:percLoad_cent
##    Data: bdta
## 
## REML criterion at convergence: 362.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.02623 -0.51313 -0.04168  0.51460  2.13594 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 13.6     3.688   
##  Residual             24.7     4.970   
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                        Estimate Std. Error t value
## (Intercept)            123.9899     1.3564   91.41
## size_pc1                 1.2045     0.5603    2.15
## Treatment.order2        -0.6514     1.3201   -0.49
## percLoad_cent           35.3986     2.2086   16.03
## percLoad2               -5.3965     9.0698   -0.59
## size_pc1:percLoad_cent  -2.3436     1.4009   -1.67
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2
## size_pc1    -0.031                            
## Trtmnt.rdr2 -0.382 -0.008                     
## percLod_cnt  0.055  0.089  0.073              
## percLoad2   -0.536  0.059 -0.183 -0.137       
## sz_pc1:prL_ -0.021  0.040 -0.181  0.031  0.280
m2.2 <- update(m2.1, .~. - size_pc1:percLoad_cent)
anova(m2.1, m2.2) ## drop all interactions
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.2: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent + 
## m2.2:     percLoad2 + (1 | Bee.ID)
## m2.1: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent + 
## m2.1:     percLoad2 + (1 | Bee.ID) + size_pc1:percLoad_cent
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## m2.2  7 395.66 410.32 -190.83   381.66                           
## m2.1  8 394.61 411.36 -189.30   378.61 3.0552      1    0.08048 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m2.2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent +  
##     percLoad2 + (1 | Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: 367.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.89123 -0.52057 -0.02436  0.57562  2.02535 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 13.86    3.723   
##  Residual             25.63    5.062   
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)      123.9447     1.3781   89.94
## size_pc1           1.2414     0.5677    2.19
## Treatment.order2  -1.0513     1.3225   -0.79
## percLoad_cent     35.5007     2.2478   15.79
## percLoad2         -1.1743     8.8665   -0.13
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2 prcLd_
## size_pc1    -0.031                     
## Trtmnt.rdr2 -0.393 -0.001              
## percLod_cnt  0.056  0.088  0.080       
## percLoad2   -0.553  0.051 -0.140 -0.152
# renaming model to simplify later typing
m2 <- m2.2

##### LRTs for main effects
## size
m3 <- update(m2, .~. - percLoad2)
anova(m2, m3, test = "Chi") # drop squared term
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m3: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent + 
## m3:     (1 | Bee.ID)
## m2: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent + 
## m2:     percLoad2 + (1 | Bee.ID)
##    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m3  6 393.68 406.24 -190.84   381.68                         
## m2  7 395.66 410.32 -190.83   381.66 0.0174      1     0.8949
summary(m3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent +  
##     (1 | Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: 373.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.87059 -0.50457 -0.03627  0.58077  2.02395 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 14.31    3.782   
##  Residual             24.73    4.973   
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)      123.8429     1.1414  108.50
## size_pc1           1.2463     0.5668    2.20
## Treatment.order2  -1.0742     1.2864   -0.84
## percLoad_cent     35.5004     2.1858   16.24
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2
## size_pc1    -0.003              
## Trtmnt.rdr2 -0.564  0.006       
## percLod_cnt -0.034  0.096  0.060
# trt order
m4 <- update(m3, .~. - Treatment.order)
anova(m3, m4) # drop treatment order
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m4: stroke.amplitude ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
## m3: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent + 
## m3:     (1 | Bee.ID)
##    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m4  5 392.42 402.89 -191.21   382.42                         
## m3  6 393.68 406.24 -190.84   381.68 0.7374      1     0.3905
summary(m4)
## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: 376.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.96941 -0.48117 -0.05909  0.48300  1.94459 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 14.53    3.811   
##  Residual             24.42    4.941   
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept)   123.3058     0.9440  130.62
## size_pc1        1.2495     0.5675    2.20
## percLoad_cent  35.6293     2.1693   16.42
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1
## size_pc1    0.000        
## percLod_cnt 0.000  0.095
# LRT for size
m5 <- update(m4, .~. - size_pc1)
anova(m4, m5) # keep size
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m5: stroke.amplitude ~ percLoad_cent + (1 | Bee.ID)
## m4: stroke.amplitude ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
##    Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)  
## m5  4 395.24 403.62 -193.62   387.24                          
## m4  5 392.42 402.89 -191.21   382.42 4.828      1      0.028 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m6 <- update(m4, .~. - percLoad_cent)
anova(m4, m6) # keep percent load
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m6: stroke.amplitude ~ size_pc1 + (1 | Bee.ID)
## m4: stroke.amplitude ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
##    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## m6  4 478.07 486.45 -235.04   470.07                             
## m4  5 392.42 402.89 -191.21   382.42 87.657      1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# summarize final model for paper
summary(m4)
## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: 376.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.96941 -0.48117 -0.05909  0.48300  1.94459 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 14.53    3.811   
##  Residual             24.42    4.941   
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept)   123.3058     0.9440  130.62
## size_pc1        1.2495     0.5675    2.20
## percLoad_cent  35.6293     2.1693   16.42
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1
## size_pc1    0.000        
## percLod_cnt 0.000  0.095
# write output
summary(m4)$coefficients  
##                 Estimate Std. Error    t value
## (Intercept)   123.305797  0.9440279 130.616690
## size_pc1        1.249463  0.5674547   2.201872
## percLoad_cent  35.629329  2.1693080  16.424283
write.csv( summary(m4)$coefficients, file = "AmpCoefs_percLoad.csv" )

strokeAmplitude diagnostics

# qq plot
qqnorm(resid(m4), main = "")
qqline(resid(m4)) # ok

# residual plot
plot(fitted(m4), resid(m3), xlab = "fitted", ylab = "residuals")
abline(0,0)

# QQPlot for group-level effects
qqnorm(ranef(m4)$Bee.ID[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(m4)$Bee.ID[[1]]) # looks good

Stroke Amplitude Visualization

preds1 <- predict(m4, newdf, re.form = NA) #predict for mean of random effect

# plot of prediction for average sized bee (line)
# raw data plotted as points
ggplot(bdta, aes(x = percLoad_cent + mean(bdta$percLoad) , y = stroke.amplitude)) + 
     geom_point(alpha = 0.3) + 
     geom_line(data = newdf, aes(x = percLoad_cent + mean(bdta$percLoad), y = preds1)) + 
     labs(x = "Mass during flight (% of starved mass in g)")

Wing Velocity

# make a full model with all two-way interactions
m1 <- lmer(wing.velocity ~  (size_pc1 +  Treatment.order + percLoad_cent)^2 + 
                percLoad2 + 
                (1|Bee.ID), data = bdta)

summary(m1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: wing.velocity ~ (size_pc1 + Treatment.order + percLoad_cent)^2 +  
##     percLoad2 + (1 | Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: -76.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.96215 -0.44762 -0.04005  0.50444  2.11529 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 0.003968 0.06299 
##  Residual             0.006025 0.07762 
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                Estimate Std. Error t value
## (Intercept)                     1.48344    0.02182   67.98
## size_pc1                       -0.04704    0.01148   -4.10
## Treatment.order2               -0.01763    0.02071   -0.85
## percLoad_cent                  -0.41560    0.06273   -6.62
## percLoad2                      -0.05676    0.14666   -0.39
## size_pc1:Treatment.order2       0.01481    0.01306    1.13
## size_pc1:percLoad_cent          0.04749    0.02300    2.06
## Treatment.order2:percLoad_cent  0.05562    0.09323    0.60
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2 s_1:T. s_1:L_
## size_pc1     0.022                                          
## Trtmnt.rdr2 -0.362  0.048                                   
## percLod_cnt  0.017  0.288  0.085                            
## percLoad2   -0.527 -0.099 -0.199 -0.153                     
## sz_pc1:Tr.2 -0.085 -0.591 -0.087 -0.287  0.240              
## sz_pc1:prL_ -0.040 -0.147 -0.200 -0.147  0.330  0.279       
## Trtmnt.2:L_  0.025 -0.242 -0.045 -0.829  0.072  0.241  0.170
### LRT's for interactions
m2.0 <- update(m1, .~. - Treatment.order:percLoad_cent )

anova(m1, m2.0)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.0: wing.velocity ~ size_pc1 + Treatment.order + percLoad_cent + 
## m2.0:     percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order + size_pc1:percLoad_cent
## m1: wing.velocity ~ (size_pc1 + Treatment.order + percLoad_cent)^2 + 
## m1:     percLoad2 + (1 | Bee.ID)
##      Df      AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## m2.0  9 -101.513 -82.664 59.756  -119.51                         
## m1   10  -99.871 -78.927 59.935  -119.87 0.3579      1     0.5497
summary(m2.0)
## Linear mixed model fit by REML ['lmerMod']
## Formula: wing.velocity ~ size_pc1 + Treatment.order + percLoad_cent +  
##     percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order + size_pc1:percLoad_cent
##    Data: bdta
## 
## REML criterion at convergence: -79.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.02111 -0.46259 -0.01217  0.49529  2.12092 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 0.004027 0.06346 
##  Residual             0.005882 0.07670 
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                           Estimate Std. Error t value
## (Intercept)                1.48322    0.02168   68.42
## size_pc1                  -0.04539    0.01109   -4.09
## Treatment.order2          -0.01704    0.02044   -0.83
## percLoad_cent             -0.38478    0.03470  -11.09
## percLoad2                 -0.06438    0.14470   -0.44
## size_pc1:Treatment.order2  0.01292    0.01253    1.03
## size_pc1:percLoad_cent     0.04506    0.02241    2.01
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2 s_1:T.
## size_pc1     0.028                                   
## Trtmnt.rdr2 -0.359  0.038                            
## percLod_cnt  0.067  0.160  0.085                     
## percLoad2   -0.528 -0.084 -0.197 -0.167              
## sz_pc1:Tr.2 -0.093 -0.561 -0.079 -0.162  0.230       
## sz_pc1:prL_ -0.045 -0.110 -0.196 -0.011  0.324  0.249
m2.1 <- update(m2.0,.~. - size_pc1:Treatment.order)
anova(m2.0, m2.1)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.1: wing.velocity ~ size_pc1 + Treatment.order + percLoad_cent + 
## m2.1:     percLoad2 + (1 | Bee.ID) + size_pc1:percLoad_cent
## m2.0: wing.velocity ~ size_pc1 + Treatment.order + percLoad_cent + 
## m2.0:     percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order + size_pc1:percLoad_cent
##      Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## m2.1  8 -102.28 -85.525 59.140  -118.28                         
## m2.0  9 -101.51 -82.664 59.756  -119.51 1.2331      1     0.2668
summary(m2.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: wing.velocity ~ size_pc1 + Treatment.order + percLoad_cent +  
##     percLoad2 + (1 | Bee.ID) + size_pc1:percLoad_cent
##    Data: bdta
## 
## REML criterion at convergence: -84.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.01405 -0.48637  0.04786  0.57874  2.09154 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 0.004062 0.06374 
##  Residual             0.005871 0.07662 
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                         Estimate Std. Error t value
## (Intercept)             1.485332   0.021598   68.77
## size_pc1               -0.038972   0.009194   -4.24
## Treatment.order2       -0.015368   0.020360   -0.75
## percLoad_cent          -0.379047   0.034213  -11.08
## percLoad2              -0.099134   0.140706   -0.70
## size_pc1:percLoad_cent  0.039276   0.021687    1.81
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2
## size_pc1    -0.029                            
## Trtmnt.rdr2 -0.368 -0.008                     
## percLod_cnt  0.052  0.084  0.073              
## percLoad2   -0.522  0.056 -0.184 -0.135       
## sz_pc1:prL_ -0.022  0.038 -0.182  0.031  0.283
m2.2 <- update(m2.1, .~. - size_pc1:percLoad_cent)
anova(m2.1, m2.2) ## drop all interactions, though this is very close
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.2: wing.velocity ~ size_pc1 + Treatment.order + percLoad_cent + 
## m2.2:     percLoad2 + (1 | Bee.ID)
## m2.1: wing.velocity ~ size_pc1 + Treatment.order + percLoad_cent + 
## m2.1:     percLoad2 + (1 | Bee.ID) + size_pc1:percLoad_cent
##      Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)  
## m2.2  7 -100.72 -86.055 57.358  -114.72                           
## m2.1  8 -102.28 -85.525 59.140  -118.28 3.5643      1    0.05903 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m2.2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: wing.velocity ~ size_pc1 + Treatment.order + percLoad_cent +  
##     percLoad2 + (1 | Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: -87.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.89330 -0.56516  0.07513  0.49055  1.93807 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 0.004108 0.06410 
##  Residual             0.006188 0.07866 
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                   Estimate Std. Error t value
## (Intercept)       1.486075   0.022031   67.45
## size_pc1         -0.039588   0.009322   -4.25
## Treatment.order2 -0.008667   0.020552   -0.42
## percLoad_cent    -0.380720   0.035078  -10.85
## percLoad2        -0.169751   0.138432   -1.23
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2 prcLd_
## size_pc1    -0.029                     
## Trtmnt.rdr2 -0.381 -0.001              
## percLod_cnt  0.054  0.084  0.080       
## percLoad2   -0.540  0.048 -0.141 -0.150
# renaming model to simplify later typing
m2 <- m2.2

##### LRTs for main effects
## size
m3 <- update(m2, .~. - Treatment.order)
anova(m2, m3, test = "Chi") # drop trt order
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m3: wing.velocity ~ size_pc1 + percLoad_cent + percLoad2 + (1 | Bee.ID)
## m2: wing.velocity ~ size_pc1 + Treatment.order + percLoad_cent + 
## m2:     percLoad2 + (1 | Bee.ID)
##    Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## m3  6 -102.52 -89.954 57.260  -114.52                         
## m2  7 -100.72 -86.055 57.358  -114.72 0.1952      1     0.6587
summary(m3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## wing.velocity ~ size_pc1 + percLoad_cent + percLoad2 + (1 | Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: -93.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.94992 -0.54108  0.04463  0.51763  1.87287 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 0.004216 0.06493 
##  Residual             0.005997 0.07744 
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                Estimate Std. Error t value
## (Intercept)    1.482722   0.020255   73.20
## size_pc1      -0.039609   0.009329   -4.25
## percLoad_cent -0.379848   0.034460  -11.02
## percLoad2     -0.179915   0.135098   -1.33
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 prcLd_
## size_pc1    -0.031              
## percLod_cnt  0.090  0.083       
## percLoad2   -0.643  0.048 -0.140
# trt order
m4 <- update(m3, .~. - percLoad2)
anova(m3, m4) # drop percLoad
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m4: wing.velocity ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
## m3: wing.velocity ~ size_pc1 + percLoad_cent + percLoad2 + (1 | Bee.ID)
##    Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## m4  5 -102.80 -92.331 56.401  -112.80                         
## m3  6 -102.52 -89.954 57.260  -114.52 1.7177      1       0.19
summary(m4)
## Linear mixed model fit by REML ['lmerMod']
## Formula: wing.velocity ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: -93.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.05997 -0.53298 -0.03019  0.59699  2.09586 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 0.003737 0.06113 
##  Residual             0.006384 0.07990 
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                Estimate Std. Error t value
## (Intercept)    1.465370   0.015198   96.42
## size_pc1      -0.038985   0.009136   -4.27
## percLoad_cent -0.385121   0.035065  -10.98
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1
## size_pc1    0.000        
## percLod_cnt 0.000  0.095
# LRT for size
m5 <- update(m4, .~. - size_pc1)
anova(m4, m5) # keep size
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m5: wing.velocity ~ percLoad_cent + (1 | Bee.ID)
## m4: wing.velocity ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
##    Df      AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)    
## m5  4  -89.649 -81.271 48.824  -97.649                             
## m4  5 -102.803 -92.331 56.401 -112.803 15.154      1  9.908e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m6 <- update(m4, .~. - percLoad_cent)
anova(m4, m6) # keep percent load
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m6: wing.velocity ~ size_pc1 + (1 | Bee.ID)
## m4: wing.velocity ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
##    Df      AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)    
## m6  4  -48.945 -40.567 28.472  -56.945                             
## m4  5 -102.803 -92.331 56.401 -112.803 55.858      1  7.791e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# summarize final model for paper
summary(m4)
## Linear mixed model fit by REML ['lmerMod']
## Formula: wing.velocity ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: -93.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.05997 -0.53298 -0.03019  0.59699  2.09586 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 0.003737 0.06113 
##  Residual             0.006384 0.07990 
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                Estimate Std. Error t value
## (Intercept)    1.465370   0.015198   96.42
## size_pc1      -0.038985   0.009136   -4.27
## percLoad_cent -0.385121   0.035065  -10.98
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1
## size_pc1    0.000        
## percLod_cnt 0.000  0.095
# write output
summary(m4)$coefficients  
##                  Estimate  Std. Error    t value
## (Intercept)    1.46537020 0.015197520  96.421669
## size_pc1      -0.03898547 0.009135556  -4.267444
## percLoad_cent -0.38512098 0.035064635 -10.983174
write.csv( summary(m4)$coefficients, file = "WingVelocity_percLoad.csv" )

Wing Velocity diagnostics

# qq plot
qqnorm(resid(m4), main = "")
qqline(resid(m4)) # good

# cook's distance
plot(influence(m4, obs = TRUE), which = 'cook')

# residual plot
plot(fitted(m4), resid(m4), xlab = "fitted", ylab = "residuals")
abline(0,0)

# QQPlot for group-level effects
qqnorm(ranef(m4)$Bee.ID[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(m4)$Bee.ID[[1]]) # looks good

Wing Velocity Visualization

preds1 <- predict(m4, newdf, re.form = NA) #predict for mean of random effect

# plot of prediction for average sized bee (line)
# raw data plotted as points
ggplot(bdta, aes(x = percLoad_cent + mean(bdta$percLoad) , y = wing.velocity)) + 
     geom_point(alpha = 0.3) + 
     geom_line(data = newdf, aes(x = percLoad_cent + mean(bdta$percLoad), y = preds1)) + 
     labs(x = "Mass during flight (% of starved mass in g)")

Part II: Which flight measurements are most closely associated with

respiratory rate in bees – wingbeat frequency, stroke amplitude, or

wing velocity?

Multilevel model approach to see the relatedness of the response variables

  • Wingbeat frequency
  • Stroke amplitude
  • Wing velocity
  • Resp. Rate
# predicting avg resp rate, including stroke amplitude, wing velocity and wingbeat frequency
# looks like we have to remove wing velocity
vif.mer(lmer(av.resp..CO2.mL.hr. ~ strokeAmp_cent + wingVel_cent + wbf_cent + percLoad_cent + Treatment.order + size_pc1 + (1|Bee.ID), data = bdta))
##   strokeAmp_cent     wingVel_cent         wbf_cent    percLoad_cent 
##       190.642035       133.592345        11.644262        10.853222 
## Treatment.order2         size_pc1 
##         1.262757         1.364694
vif.mer(lmer(av.resp..CO2.mL.hr. ~ strokeAmp_cent  + wbf_cent + percLoad_cent + Treatment.order + size_pc1 + (1|Bee.ID), data = bdta)) # still worryingly high for percLoad and stroke amplitude
##   strokeAmp_cent         wbf_cent    percLoad_cent Treatment.order2 
##         7.741480         1.830209         8.182520         1.255652 
##         size_pc1 
##         1.355159
mm1 <- lmer(av.resp..CO2.mL.hr. ~ (strokeAmp_cent + wbf_cent + 
                                        percLoad_cent + Treatment.order + 
                                        size_pc1)^2 + 
                 strokeAmp2 + wbf2 + percLoad2 + 
                 (1|Bee.ID), data = bdta)
## Warning: Some predictor variables are on very different scales: consider
## rescaling
summary(mm1)  # warning says to rescale
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## av.resp..CO2.mL.hr. ~ (strokeAmp_cent + wbf_cent + percLoad_cent +  
##     Treatment.order + size_pc1)^2 + strokeAmp2 + wbf2 + percLoad2 +  
##     (1 | Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: 218.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.55796 -0.28776 -0.05692  0.20517  2.08071 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 1.1866   1.0893  
##  Residual             0.2396   0.4895  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                  Estimate Std. Error t value
## (Intercept)                     10.265659   0.344190  29.826
## strokeAmp_cent                  -0.116525   0.032532  -3.582
## wbf_cent                         0.160723   0.036796   4.368
## percLoad_cent                    8.410281   1.344833   6.254
## Treatment.order2                 0.052521   0.163024   0.322
## size_pc1                         1.936582   0.189350  10.228
## strokeAmp2                       0.007777   0.002598   2.994
## wbf2                            -0.002007   0.003538  -0.567
## percLoad2                       12.861989   4.588437   2.803
## strokeAmp_cent:wbf_cent          0.005223   0.004034   1.295
## strokeAmp_cent:percLoad_cent    -0.712136   0.201023  -3.543
## strokeAmp_cent:Treatment.order2  0.093952   0.038242   2.457
## strokeAmp_cent:size_pc1          0.020108   0.014717   1.366
## wbf_cent:percLoad_cent          -0.225045   0.186731  -1.205
## wbf_cent:Treatment.order2       -0.070432   0.036113  -1.950
## wbf_cent:size_pc1               -0.042046   0.018814  -2.235
## percLoad_cent:Treatment.order2  -3.626144   1.421979  -2.550
## percLoad_cent:size_pc1          -0.467014   0.594009  -0.786
## Treatment.order2:size_pc1       -0.575909   0.150526  -3.826
## 
## Correlation matrix not shown by default, as p = 19 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
# rescale numeric predictors
bdta_scaled = scale(bdta[, c("av.resp..CO2.mL.hr.", 
                             "strokeAmp_cent", 
                             "wbf_cent",
                             "percLoad_cent",
                             "size_pc1",
                             "strokeAmp2",
                             "wbf2",
                             "percLoad2")], 
                    center = FALSE, 
                    scale = TRUE)


colnames(bdta_scaled) <- paste0(colnames(bdta_scaled), "_scaled")
scatterplotMatrix(bdta_scaled)

bdta = cbind(bdta, bdta_scaled)

mm1 <- lmer(av.resp..CO2.mL.hr. ~ (strokeAmp_cent_scaled + wbf_cent_scaled + 
                                        percLoad_cent_scaled + Treatment.order + 
                                        size_pc1_scaled)^2 + 
                 strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + 
                 (1|Bee.ID), data = bdta)
summary(mm1)  # warning says to rescale
## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ (strokeAmp_cent_scaled + wbf_cent_scaled +  
##     percLoad_cent_scaled + Treatment.order + size_pc1_scaled)^2 +  
##     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |      Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: 164.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.55796 -0.28776 -0.05692  0.20517  2.08071 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 1.1866   1.0893  
##  Residual             0.2396   0.4895  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                            Estimate Std. Error t value
## (Intercept)                                10.26566    0.34419  29.826
## strokeAmp_cent_scaled                      -1.43098    0.39951  -3.582
## wbf_cent_scaled                             1.22494    0.28044   4.368
## percLoad_cent_scaled                        2.63390    0.42117   6.254
## Treatment.order2                            0.05252    0.16302   0.322
## size_pc1_scaled                             3.26359    0.31910  10.228
## strokeAmp2_scaled                           1.68151    0.56164   2.994
## wbf2_scaled                                -0.16967    0.29913  -0.567
## percLoad2_scaled                            1.61846    0.57738   2.803
## strokeAmp_cent_scaled:wbf_cent_scaled       0.48889    0.37760   1.295
## strokeAmp_cent_scaled:percLoad_cent_scaled -2.73884    0.77312  -3.543
## strokeAmp_cent_scaled:Treatment.order2      1.15377    0.46962   2.457
## strokeAmp_cent_scaled:size_pc1_scaled       0.41615    0.30457   1.366
## wbf_cent_scaled:percLoad_cent_scaled       -0.53715    0.44570  -1.205
## wbf_cent_scaled:Treatment.order2           -0.53679    0.27524  -1.950
## wbf_cent_scaled:size_pc1_scaled            -0.54003    0.24164  -2.235
## percLoad_cent_scaled:Treatment.order2      -1.13562    0.44533  -2.550
## percLoad_cent_scaled:size_pc1_scaled       -0.24648    0.31350  -0.786
## Treatment.order2:size_pc1_scaled           -0.97054    0.25367  -3.826
## 
## Correlation matrix not shown by default, as p = 19 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it
# LRT's
mm2 <- update(mm1, .~. -  percLoad_cent_scaled:size_pc1_scaled)
anova(mm1, mm2)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## mm2: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled + 
## mm2:     percLoad_cent_scaled + Treatment.order + size_pc1_scaled + 
## mm2:     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 | 
## mm2:     Bee.ID) + strokeAmp_cent_scaled:wbf_cent_scaled + strokeAmp_cent_scaled:percLoad_cent_scaled + 
## mm2:     strokeAmp_cent_scaled:Treatment.order + strokeAmp_cent_scaled:size_pc1_scaled + 
## mm2:     wbf_cent_scaled:percLoad_cent_scaled + wbf_cent_scaled:Treatment.order + 
## mm2:     wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order + 
## mm2:     Treatment.order:size_pc1_scaled
## mm1: av.resp..CO2.mL.hr. ~ (strokeAmp_cent_scaled + wbf_cent_scaled + 
## mm1:     percLoad_cent_scaled + Treatment.order + size_pc1_scaled)^2 + 
## mm1:     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 | 
## mm1:     Bee.ID)
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## mm2 20 172.68 214.56 -66.338   132.68                         
## mm1 21 172.85 216.83 -65.425   130.85 1.8248      1     0.1767
summary(mm2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +  
##     percLoad_cent_scaled + Treatment.order + size_pc1_scaled +  
##     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |  
##     Bee.ID) + strokeAmp_cent_scaled:wbf_cent_scaled + strokeAmp_cent_scaled:percLoad_cent_scaled +  
##     strokeAmp_cent_scaled:Treatment.order + strokeAmp_cent_scaled:size_pc1_scaled +  
##     wbf_cent_scaled:percLoad_cent_scaled + wbf_cent_scaled:Treatment.order +  
##     wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order +  
##     Treatment.order:size_pc1_scaled
##    Data: bdta
## 
## REML criterion at convergence: 164.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.64217 -0.33793 -0.05746  0.20194  2.17217 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 1.1457   1.0704  
##  Residual             0.2449   0.4949  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                            Estimate Std. Error t value
## (Intercept)                                10.26035    0.34280  29.931
## strokeAmp_cent_scaled                      -1.38524    0.39805  -3.480
## wbf_cent_scaled                             1.25452    0.27883   4.499
## percLoad_cent_scaled                        2.58315    0.41822   6.177
## Treatment.order2                            0.05135    0.16419   0.313
## size_pc1_scaled                             3.29213    0.31545  10.436
## strokeAmp2_scaled                           1.62369    0.56341   2.882
## wbf2_scaled                                -0.21584    0.29498  -0.732
## percLoad2_scaled                            1.50825    0.56769   2.657
## strokeAmp_cent_scaled:wbf_cent_scaled       0.35059    0.33972   1.032
## strokeAmp_cent_scaled:percLoad_cent_scaled -2.57308    0.75769  -3.396
## strokeAmp_cent_scaled:Treatment.order2      1.10875    0.47154   2.351
## strokeAmp_cent_scaled:size_pc1_scaled       0.20581    0.15000   1.372
## wbf_cent_scaled:percLoad_cent_scaled       -0.33822    0.37227  -0.909
## wbf_cent_scaled:Treatment.order2           -0.52811    0.27740  -1.904
## wbf_cent_scaled:size_pc1_scaled            -0.57156    0.23678  -2.414
## percLoad_cent_scaled:Treatment.order2      -1.09799    0.44528  -2.466
## Treatment.order2:size_pc1_scaled           -0.98374    0.25484  -3.860
## 
## Correlation matrix not shown by default, as p = 18 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it
mm3 <- update(mm2, .~. - wbf_cent_scaled:percLoad_cent_scaled)
anova(mm2,mm3 )
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## mm3: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled + 
## mm3:     percLoad_cent_scaled + Treatment.order + size_pc1_scaled + 
## mm3:     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 | 
## mm3:     Bee.ID) + strokeAmp_cent_scaled:wbf_cent_scaled + strokeAmp_cent_scaled:percLoad_cent_scaled + 
## mm3:     strokeAmp_cent_scaled:Treatment.order + strokeAmp_cent_scaled:size_pc1_scaled + 
## mm3:     wbf_cent_scaled:Treatment.order + wbf_cent_scaled:size_pc1_scaled + 
## mm3:     percLoad_cent_scaled:Treatment.order + Treatment.order:size_pc1_scaled
## mm2: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled + 
## mm2:     percLoad_cent_scaled + Treatment.order + size_pc1_scaled + 
## mm2:     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 | 
## mm2:     Bee.ID) + strokeAmp_cent_scaled:wbf_cent_scaled + strokeAmp_cent_scaled:percLoad_cent_scaled + 
## mm2:     strokeAmp_cent_scaled:Treatment.order + strokeAmp_cent_scaled:size_pc1_scaled + 
## mm2:     wbf_cent_scaled:percLoad_cent_scaled + wbf_cent_scaled:Treatment.order + 
## mm2:     wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order + 
## mm2:     Treatment.order:size_pc1_scaled
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## mm3 19 172.39 212.18 -67.193   134.39                         
## mm2 20 172.68 214.56 -66.338   132.68 1.7113      1     0.1908
summary(mm3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +  
##     percLoad_cent_scaled + Treatment.order + size_pc1_scaled +  
##     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |  
##     Bee.ID) + strokeAmp_cent_scaled:wbf_cent_scaled + strokeAmp_cent_scaled:percLoad_cent_scaled +  
##     strokeAmp_cent_scaled:Treatment.order + strokeAmp_cent_scaled:size_pc1_scaled +  
##     wbf_cent_scaled:Treatment.order + wbf_cent_scaled:size_pc1_scaled +  
##     percLoad_cent_scaled:Treatment.order + Treatment.order:size_pc1_scaled
##    Data: bdta
## 
## REML criterion at convergence: 165
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.64387 -0.33501 -0.04398  0.30860  2.16600 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 1.1202   1.0584  
##  Residual             0.2491   0.4991  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                            Estimate Std. Error t value
## (Intercept)                                10.32638    0.33546  30.783
## strokeAmp_cent_scaled                      -1.24881    0.37193  -3.358
## wbf_cent_scaled                             1.24684    0.27979   4.456
## percLoad_cent_scaled                        2.45007    0.39328   6.230
## Treatment.order2                            0.08464    0.16095   0.526
## size_pc1_scaled                             3.22835    0.30701  10.515
## strokeAmp2_scaled                           1.37824    0.50568   2.726
## wbf2_scaled                                -0.32636    0.27002  -1.209
## percLoad2_scaled                            1.14468    0.41345   2.769
## strokeAmp_cent_scaled:wbf_cent_scaled       0.06994    0.14601   0.479
## strokeAmp_cent_scaled:percLoad_cent_scaled -2.12997    0.59831  -3.560
## strokeAmp_cent_scaled:Treatment.order2      1.00405    0.46268   2.170
## strokeAmp_cent_scaled:size_pc1_scaled       0.24015    0.14555   1.650
## wbf_cent_scaled:Treatment.order2           -0.52241    0.27921  -1.871
## wbf_cent_scaled:size_pc1_scaled            -0.56075    0.23620  -2.374
## percLoad_cent_scaled:Treatment.order2      -0.96560    0.42207  -2.288
## Treatment.order2:size_pc1_scaled           -0.92574    0.24880  -3.721
## 
## Correlation matrix not shown by default, as p = 17 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it
mm4 <- update(mm3, .~.  - strokeAmp_cent_scaled:size_pc1_scaled)
anova(mm3, mm4) # keep interaction
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## mm4: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled + 
## mm4:     percLoad_cent_scaled + Treatment.order + size_pc1_scaled + 
## mm4:     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 | 
## mm4:     Bee.ID) + strokeAmp_cent_scaled:wbf_cent_scaled + strokeAmp_cent_scaled:percLoad_cent_scaled + 
## mm4:     strokeAmp_cent_scaled:Treatment.order + wbf_cent_scaled:Treatment.order + 
## mm4:     wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order + 
## mm4:     Treatment.order:size_pc1_scaled
## mm3: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled + 
## mm3:     percLoad_cent_scaled + Treatment.order + size_pc1_scaled + 
## mm3:     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 | 
## mm3:     Bee.ID) + strokeAmp_cent_scaled:wbf_cent_scaled + strokeAmp_cent_scaled:percLoad_cent_scaled + 
## mm3:     strokeAmp_cent_scaled:Treatment.order + strokeAmp_cent_scaled:size_pc1_scaled + 
## mm3:     wbf_cent_scaled:Treatment.order + wbf_cent_scaled:size_pc1_scaled + 
## mm3:     percLoad_cent_scaled:Treatment.order + Treatment.order:size_pc1_scaled
##     Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)  
## mm4 18 175.30 213.00 -69.651   139.30                          
## mm3 19 172.39 212.18 -67.193   134.39 4.916      1    0.02661 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mm3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +  
##     percLoad_cent_scaled + Treatment.order + size_pc1_scaled +  
##     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |  
##     Bee.ID) + strokeAmp_cent_scaled:wbf_cent_scaled + strokeAmp_cent_scaled:percLoad_cent_scaled +  
##     strokeAmp_cent_scaled:Treatment.order + strokeAmp_cent_scaled:size_pc1_scaled +  
##     wbf_cent_scaled:Treatment.order + wbf_cent_scaled:size_pc1_scaled +  
##     percLoad_cent_scaled:Treatment.order + Treatment.order:size_pc1_scaled
##    Data: bdta
## 
## REML criterion at convergence: 165
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.64387 -0.33501 -0.04398  0.30860  2.16600 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 1.1202   1.0584  
##  Residual             0.2491   0.4991  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                            Estimate Std. Error t value
## (Intercept)                                10.32638    0.33546  30.783
## strokeAmp_cent_scaled                      -1.24881    0.37193  -3.358
## wbf_cent_scaled                             1.24684    0.27979   4.456
## percLoad_cent_scaled                        2.45007    0.39328   6.230
## Treatment.order2                            0.08464    0.16095   0.526
## size_pc1_scaled                             3.22835    0.30701  10.515
## strokeAmp2_scaled                           1.37824    0.50568   2.726
## wbf2_scaled                                -0.32636    0.27002  -1.209
## percLoad2_scaled                            1.14468    0.41345   2.769
## strokeAmp_cent_scaled:wbf_cent_scaled       0.06994    0.14601   0.479
## strokeAmp_cent_scaled:percLoad_cent_scaled -2.12997    0.59831  -3.560
## strokeAmp_cent_scaled:Treatment.order2      1.00405    0.46268   2.170
## strokeAmp_cent_scaled:size_pc1_scaled       0.24015    0.14555   1.650
## wbf_cent_scaled:Treatment.order2           -0.52241    0.27921  -1.871
## wbf_cent_scaled:size_pc1_scaled            -0.56075    0.23620  -2.374
## percLoad_cent_scaled:Treatment.order2      -0.96560    0.42207  -2.288
## Treatment.order2:size_pc1_scaled           -0.92574    0.24880  -3.721
## 
## Correlation matrix not shown by default, as p = 17 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it
mm5 <- update(mm3, .~. - wbf_cent_scaled:Treatment.order)
anova(mm3, mm5) # keep interaction
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## mm5: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled + 
## mm5:     percLoad_cent_scaled + Treatment.order + size_pc1_scaled + 
## mm5:     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 | 
## mm5:     Bee.ID) + strokeAmp_cent_scaled:wbf_cent_scaled + strokeAmp_cent_scaled:percLoad_cent_scaled + 
## mm5:     strokeAmp_cent_scaled:Treatment.order + strokeAmp_cent_scaled:size_pc1_scaled + 
## mm5:     wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order + 
## mm5:     Treatment.order:size_pc1_scaled
## mm3: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled + 
## mm3:     percLoad_cent_scaled + Treatment.order + size_pc1_scaled + 
## mm3:     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 | 
## mm3:     Bee.ID) + strokeAmp_cent_scaled:wbf_cent_scaled + strokeAmp_cent_scaled:percLoad_cent_scaled + 
## mm3:     strokeAmp_cent_scaled:Treatment.order + strokeAmp_cent_scaled:size_pc1_scaled + 
## mm3:     wbf_cent_scaled:Treatment.order + wbf_cent_scaled:size_pc1_scaled + 
## mm3:     percLoad_cent_scaled:Treatment.order + Treatment.order:size_pc1_scaled
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## mm5 18 176.32 214.02 -70.159   140.32                           
## mm3 19 172.39 212.18 -67.193   134.39 5.9309      1    0.01488 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mm3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +  
##     percLoad_cent_scaled + Treatment.order + size_pc1_scaled +  
##     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |  
##     Bee.ID) + strokeAmp_cent_scaled:wbf_cent_scaled + strokeAmp_cent_scaled:percLoad_cent_scaled +  
##     strokeAmp_cent_scaled:Treatment.order + strokeAmp_cent_scaled:size_pc1_scaled +  
##     wbf_cent_scaled:Treatment.order + wbf_cent_scaled:size_pc1_scaled +  
##     percLoad_cent_scaled:Treatment.order + Treatment.order:size_pc1_scaled
##    Data: bdta
## 
## REML criterion at convergence: 165
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.64387 -0.33501 -0.04398  0.30860  2.16600 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 1.1202   1.0584  
##  Residual             0.2491   0.4991  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                            Estimate Std. Error t value
## (Intercept)                                10.32638    0.33546  30.783
## strokeAmp_cent_scaled                      -1.24881    0.37193  -3.358
## wbf_cent_scaled                             1.24684    0.27979   4.456
## percLoad_cent_scaled                        2.45007    0.39328   6.230
## Treatment.order2                            0.08464    0.16095   0.526
## size_pc1_scaled                             3.22835    0.30701  10.515
## strokeAmp2_scaled                           1.37824    0.50568   2.726
## wbf2_scaled                                -0.32636    0.27002  -1.209
## percLoad2_scaled                            1.14468    0.41345   2.769
## strokeAmp_cent_scaled:wbf_cent_scaled       0.06994    0.14601   0.479
## strokeAmp_cent_scaled:percLoad_cent_scaled -2.12997    0.59831  -3.560
## strokeAmp_cent_scaled:Treatment.order2      1.00405    0.46268   2.170
## strokeAmp_cent_scaled:size_pc1_scaled       0.24015    0.14555   1.650
## wbf_cent_scaled:Treatment.order2           -0.52241    0.27921  -1.871
## wbf_cent_scaled:size_pc1_scaled            -0.56075    0.23620  -2.374
## percLoad_cent_scaled:Treatment.order2      -0.96560    0.42207  -2.288
## Treatment.order2:size_pc1_scaled           -0.92574    0.24880  -3.721
## 
## Correlation matrix not shown by default, as p = 17 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it
mm6 <- update(mm3, .~. - strokeAmp_cent_scaled:wbf_cent_scaled)
anova(mm3, mm6) # drop interaction
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## mm6: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled + 
## mm6:     percLoad_cent_scaled + Treatment.order + size_pc1_scaled + 
## mm6:     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 | 
## mm6:     Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order + 
## mm6:     strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order + 
## mm6:     wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order + 
## mm6:     Treatment.order:size_pc1_scaled
## mm3: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled + 
## mm3:     percLoad_cent_scaled + Treatment.order + size_pc1_scaled + 
## mm3:     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 | 
## mm3:     Bee.ID) + strokeAmp_cent_scaled:wbf_cent_scaled + strokeAmp_cent_scaled:percLoad_cent_scaled + 
## mm3:     strokeAmp_cent_scaled:Treatment.order + strokeAmp_cent_scaled:size_pc1_scaled + 
## mm3:     wbf_cent_scaled:Treatment.order + wbf_cent_scaled:size_pc1_scaled + 
## mm3:     percLoad_cent_scaled:Treatment.order + Treatment.order:size_pc1_scaled
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## mm6 18 171.04 208.74 -67.518   135.04                         
## mm3 19 172.39 212.18 -67.193   134.39 0.6498      1     0.4202
summary(mm6)
## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +  
##     percLoad_cent_scaled + Treatment.order + size_pc1_scaled +  
##     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |  
##     Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order +  
##     strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order +  
##     wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order +  
##     Treatment.order:size_pc1_scaled
##    Data: bdta
## 
## REML criterion at convergence: 163.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6717 -0.3034 -0.0680  0.3528  2.2000 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 1.1046   1.0510  
##  Residual             0.2438   0.4937  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                            Estimate Std. Error t value
## (Intercept)                                10.33167    0.33224  31.097
## strokeAmp_cent_scaled                      -1.21537    0.36143  -3.363
## wbf_cent_scaled                             1.20177    0.26092   4.606
## percLoad_cent_scaled                        2.45005    0.38932   6.293
## Treatment.order2                            0.06881    0.15588   0.441
## size_pc1_scaled                             3.18673    0.29199  10.914
## strokeAmp2_scaled                           1.36399    0.49926   2.732
## wbf2_scaled                                -0.27572    0.24625  -1.120
## percLoad2_scaled                            1.12004    0.40566   2.761
## strokeAmp_cent_scaled:percLoad_cent_scaled -2.08682    0.58439  -3.571
## strokeAmp_cent_scaled:Treatment.order2      0.96049    0.44821   2.143
## strokeAmp_cent_scaled:size_pc1_scaled       0.19023    0.10033   1.896
## wbf_cent_scaled:Treatment.order2           -0.48672    0.26616  -1.829
## wbf_cent_scaled:size_pc1_scaled            -0.51630    0.21509  -2.400
## percLoad_cent_scaled:Treatment.order2      -0.98064    0.41651  -2.354
## Treatment.order2:size_pc1_scaled           -0.88736    0.23299  -3.809
## 
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it
mm7 <- update(mm6, .~. - wbf_cent_scaled:size_pc1_scaled)
anova(mm6, mm7) # keep
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## mm7: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled + 
## mm7:     percLoad_cent_scaled + Treatment.order + size_pc1_scaled + 
## mm7:     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 | 
## mm7:     Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order + 
## mm7:     strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order + 
## mm7:     percLoad_cent_scaled:Treatment.order + Treatment.order:size_pc1_scaled
## mm6: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled + 
## mm6:     percLoad_cent_scaled + Treatment.order + size_pc1_scaled + 
## mm6:     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 | 
## mm6:     Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order + 
## mm6:     strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order + 
## mm6:     wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order + 
## mm6:     Treatment.order:size_pc1_scaled
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)   
## mm7 17 176.71 212.32 -71.357   142.71                            
## mm6 18 171.04 208.74 -67.518   135.04 7.6766      1   0.005594 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mm6)
## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +  
##     percLoad_cent_scaled + Treatment.order + size_pc1_scaled +  
##     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |  
##     Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order +  
##     strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order +  
##     wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order +  
##     Treatment.order:size_pc1_scaled
##    Data: bdta
## 
## REML criterion at convergence: 163.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6717 -0.3034 -0.0680  0.3528  2.2000 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 1.1046   1.0510  
##  Residual             0.2438   0.4937  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                            Estimate Std. Error t value
## (Intercept)                                10.33167    0.33224  31.097
## strokeAmp_cent_scaled                      -1.21537    0.36143  -3.363
## wbf_cent_scaled                             1.20177    0.26092   4.606
## percLoad_cent_scaled                        2.45005    0.38932   6.293
## Treatment.order2                            0.06881    0.15588   0.441
## size_pc1_scaled                             3.18673    0.29199  10.914
## strokeAmp2_scaled                           1.36399    0.49926   2.732
## wbf2_scaled                                -0.27572    0.24625  -1.120
## percLoad2_scaled                            1.12004    0.40566   2.761
## strokeAmp_cent_scaled:percLoad_cent_scaled -2.08682    0.58439  -3.571
## strokeAmp_cent_scaled:Treatment.order2      0.96049    0.44821   2.143
## strokeAmp_cent_scaled:size_pc1_scaled       0.19023    0.10033   1.896
## wbf_cent_scaled:Treatment.order2           -0.48672    0.26616  -1.829
## wbf_cent_scaled:size_pc1_scaled            -0.51630    0.21509  -2.400
## percLoad_cent_scaled:Treatment.order2      -0.98064    0.41651  -2.354
## Treatment.order2:size_pc1_scaled           -0.88736    0.23299  -3.809
## 
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it
mm8 <- update(mm6, .~.  - percLoad_cent_scaled:Treatment.order)
anova(mm6, mm8) # keep
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## mm8: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled + 
## mm8:     percLoad_cent_scaled + Treatment.order + size_pc1_scaled + 
## mm8:     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 | 
## mm8:     Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order + 
## mm8:     strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order + 
## mm8:     wbf_cent_scaled:size_pc1_scaled + Treatment.order:size_pc1_scaled
## mm6: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled + 
## mm6:     percLoad_cent_scaled + Treatment.order + size_pc1_scaled + 
## mm6:     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 | 
## mm6:     Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order + 
## mm6:     strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order + 
## mm6:     wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order + 
## mm6:     Treatment.order:size_pc1_scaled
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)   
## mm8 17 176.86 212.47 -71.431   142.86                            
## mm6 18 171.04 208.74 -67.518   135.04 7.8259      1    0.00515 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mm6)
## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +  
##     percLoad_cent_scaled + Treatment.order + size_pc1_scaled +  
##     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |  
##     Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order +  
##     strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order +  
##     wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order +  
##     Treatment.order:size_pc1_scaled
##    Data: bdta
## 
## REML criterion at convergence: 163.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6717 -0.3034 -0.0680  0.3528  2.2000 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 1.1046   1.0510  
##  Residual             0.2438   0.4937  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                            Estimate Std. Error t value
## (Intercept)                                10.33167    0.33224  31.097
## strokeAmp_cent_scaled                      -1.21537    0.36143  -3.363
## wbf_cent_scaled                             1.20177    0.26092   4.606
## percLoad_cent_scaled                        2.45005    0.38932   6.293
## Treatment.order2                            0.06881    0.15588   0.441
## size_pc1_scaled                             3.18673    0.29199  10.914
## strokeAmp2_scaled                           1.36399    0.49926   2.732
## wbf2_scaled                                -0.27572    0.24625  -1.120
## percLoad2_scaled                            1.12004    0.40566   2.761
## strokeAmp_cent_scaled:percLoad_cent_scaled -2.08682    0.58439  -3.571
## strokeAmp_cent_scaled:Treatment.order2      0.96049    0.44821   2.143
## strokeAmp_cent_scaled:size_pc1_scaled       0.19023    0.10033   1.896
## wbf_cent_scaled:Treatment.order2           -0.48672    0.26616  -1.829
## wbf_cent_scaled:size_pc1_scaled            -0.51630    0.21509  -2.400
## percLoad_cent_scaled:Treatment.order2      -0.98064    0.41651  -2.354
## Treatment.order2:size_pc1_scaled           -0.88736    0.23299  -3.809
## 
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it
mm9 <- update(mm6, .~. - Treatment.order:size_pc1_scaled)
anova(mm6, mm9) # keep
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## mm9: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled + 
## mm9:     percLoad_cent_scaled + Treatment.order + size_pc1_scaled + 
## mm9:     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 | 
## mm9:     Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order + 
## mm9:     strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order + 
## mm9:     wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order
## mm6: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled + 
## mm6:     percLoad_cent_scaled + Treatment.order + size_pc1_scaled + 
## mm6:     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 | 
## mm6:     Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order + 
## mm6:     strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order + 
## mm6:     wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order + 
## mm6:     Treatment.order:size_pc1_scaled
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## mm9 17 186.96 222.57 -76.481   152.96                             
## mm6 18 171.04 208.74 -67.518   135.04 17.926      1  2.296e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mm6)
## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +  
##     percLoad_cent_scaled + Treatment.order + size_pc1_scaled +  
##     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |  
##     Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order +  
##     strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order +  
##     wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order +  
##     Treatment.order:size_pc1_scaled
##    Data: bdta
## 
## REML criterion at convergence: 163.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6717 -0.3034 -0.0680  0.3528  2.2000 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 1.1046   1.0510  
##  Residual             0.2438   0.4937  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                            Estimate Std. Error t value
## (Intercept)                                10.33167    0.33224  31.097
## strokeAmp_cent_scaled                      -1.21537    0.36143  -3.363
## wbf_cent_scaled                             1.20177    0.26092   4.606
## percLoad_cent_scaled                        2.45005    0.38932   6.293
## Treatment.order2                            0.06881    0.15588   0.441
## size_pc1_scaled                             3.18673    0.29199  10.914
## strokeAmp2_scaled                           1.36399    0.49926   2.732
## wbf2_scaled                                -0.27572    0.24625  -1.120
## percLoad2_scaled                            1.12004    0.40566   2.761
## strokeAmp_cent_scaled:percLoad_cent_scaled -2.08682    0.58439  -3.571
## strokeAmp_cent_scaled:Treatment.order2      0.96049    0.44821   2.143
## strokeAmp_cent_scaled:size_pc1_scaled       0.19023    0.10033   1.896
## wbf_cent_scaled:Treatment.order2           -0.48672    0.26616  -1.829
## wbf_cent_scaled:size_pc1_scaled            -0.51630    0.21509  -2.400
## percLoad_cent_scaled:Treatment.order2      -0.98064    0.41651  -2.354
## Treatment.order2:size_pc1_scaled           -0.88736    0.23299  -3.809
## 
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it
mm10 <- update(mm6 , .~. - strokeAmp_cent_scaled:Treatment.order)
anova(mm6, mm10) # keep interaction
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## mm10: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled + 
## mm10:     percLoad_cent_scaled + Treatment.order + size_pc1_scaled + 
## mm10:     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 | 
## mm10:     Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:size_pc1_scaled + 
## mm10:     wbf_cent_scaled:Treatment.order + wbf_cent_scaled:size_pc1_scaled + 
## mm10:     percLoad_cent_scaled:Treatment.order + Treatment.order:size_pc1_scaled
## mm6: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled + 
## mm6:     percLoad_cent_scaled + Treatment.order + size_pc1_scaled + 
## mm6:     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 | 
## mm6:     Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order + 
## mm6:     strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order + 
## mm6:     wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order + 
## mm6:     Treatment.order:size_pc1_scaled
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)   
## mm10 17 176.19 211.79 -71.093   142.19                            
## mm6  18 171.04 208.74 -67.518   135.04 7.1495      1   0.007499 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mm11 <- update(mm6, .~. - strokeAmp_cent_scaled:Treatment.order)
anova(mm6, mm11) # keep interaction
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## mm11: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled + 
## mm11:     percLoad_cent_scaled + Treatment.order + size_pc1_scaled + 
## mm11:     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 | 
## mm11:     Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:size_pc1_scaled + 
## mm11:     wbf_cent_scaled:Treatment.order + wbf_cent_scaled:size_pc1_scaled + 
## mm11:     percLoad_cent_scaled:Treatment.order + Treatment.order:size_pc1_scaled
## mm6: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled + 
## mm6:     percLoad_cent_scaled + Treatment.order + size_pc1_scaled + 
## mm6:     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 | 
## mm6:     Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order + 
## mm6:     strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order + 
## mm6:     wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order + 
## mm6:     Treatment.order:size_pc1_scaled
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)   
## mm11 17 176.19 211.79 -71.093   142.19                            
## mm6  18 171.04 208.74 -67.518   135.04 7.1495      1   0.007499 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mm6)
## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +  
##     percLoad_cent_scaled + Treatment.order + size_pc1_scaled +  
##     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |  
##     Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order +  
##     strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order +  
##     wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order +  
##     Treatment.order:size_pc1_scaled
##    Data: bdta
## 
## REML criterion at convergence: 163.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6717 -0.3034 -0.0680  0.3528  2.2000 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 1.1046   1.0510  
##  Residual             0.2438   0.4937  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                            Estimate Std. Error t value
## (Intercept)                                10.33167    0.33224  31.097
## strokeAmp_cent_scaled                      -1.21537    0.36143  -3.363
## wbf_cent_scaled                             1.20177    0.26092   4.606
## percLoad_cent_scaled                        2.45005    0.38932   6.293
## Treatment.order2                            0.06881    0.15588   0.441
## size_pc1_scaled                             3.18673    0.29199  10.914
## strokeAmp2_scaled                           1.36399    0.49926   2.732
## wbf2_scaled                                -0.27572    0.24625  -1.120
## percLoad2_scaled                            1.12004    0.40566   2.761
## strokeAmp_cent_scaled:percLoad_cent_scaled -2.08682    0.58439  -3.571
## strokeAmp_cent_scaled:Treatment.order2      0.96049    0.44821   2.143
## strokeAmp_cent_scaled:size_pc1_scaled       0.19023    0.10033   1.896
## wbf_cent_scaled:Treatment.order2           -0.48672    0.26616  -1.829
## wbf_cent_scaled:size_pc1_scaled            -0.51630    0.21509  -2.400
## percLoad_cent_scaled:Treatment.order2      -0.98064    0.41651  -2.354
## Treatment.order2:size_pc1_scaled           -0.88736    0.23299  -3.809
## 
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it
# move on to squared terms
mm12 <- update(mm6 , .~. - wbf2_scaled)
anova(mm6, mm11) # keep wbf2
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## mm11: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled + 
## mm11:     percLoad_cent_scaled + Treatment.order + size_pc1_scaled + 
## mm11:     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 | 
## mm11:     Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:size_pc1_scaled + 
## mm11:     wbf_cent_scaled:Treatment.order + wbf_cent_scaled:size_pc1_scaled + 
## mm11:     percLoad_cent_scaled:Treatment.order + Treatment.order:size_pc1_scaled
## mm6: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled + 
## mm6:     percLoad_cent_scaled + Treatment.order + size_pc1_scaled + 
## mm6:     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 | 
## mm6:     Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order + 
## mm6:     strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order + 
## mm6:     wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order + 
## mm6:     Treatment.order:size_pc1_scaled
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)   
## mm11 17 176.19 211.79 -71.093   142.19                            
## mm6  18 171.04 208.74 -67.518   135.04 7.1495      1   0.007499 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mm6)
## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +  
##     percLoad_cent_scaled + Treatment.order + size_pc1_scaled +  
##     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |  
##     Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order +  
##     strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order +  
##     wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order +  
##     Treatment.order:size_pc1_scaled
##    Data: bdta
## 
## REML criterion at convergence: 163.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6717 -0.3034 -0.0680  0.3528  2.2000 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 1.1046   1.0510  
##  Residual             0.2438   0.4937  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                            Estimate Std. Error t value
## (Intercept)                                10.33167    0.33224  31.097
## strokeAmp_cent_scaled                      -1.21537    0.36143  -3.363
## wbf_cent_scaled                             1.20177    0.26092   4.606
## percLoad_cent_scaled                        2.45005    0.38932   6.293
## Treatment.order2                            0.06881    0.15588   0.441
## size_pc1_scaled                             3.18673    0.29199  10.914
## strokeAmp2_scaled                           1.36399    0.49926   2.732
## wbf2_scaled                                -0.27572    0.24625  -1.120
## percLoad2_scaled                            1.12004    0.40566   2.761
## strokeAmp_cent_scaled:percLoad_cent_scaled -2.08682    0.58439  -3.571
## strokeAmp_cent_scaled:Treatment.order2      0.96049    0.44821   2.143
## strokeAmp_cent_scaled:size_pc1_scaled       0.19023    0.10033   1.896
## wbf_cent_scaled:Treatment.order2           -0.48672    0.26616  -1.829
## wbf_cent_scaled:size_pc1_scaled            -0.51630    0.21509  -2.400
## percLoad_cent_scaled:Treatment.order2      -0.98064    0.41651  -2.354
## Treatment.order2:size_pc1_scaled           -0.88736    0.23299  -3.809
## 
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it
mm13 <- update(mm6, .~. - percLoad2_scaled)
anova(mm6, mm13) # keep
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## mm13: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled + 
## mm13:     percLoad_cent_scaled + Treatment.order + size_pc1_scaled + 
## mm13:     strokeAmp2_scaled + wbf2_scaled + (1 | Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + 
## mm13:     strokeAmp_cent_scaled:Treatment.order + strokeAmp_cent_scaled:size_pc1_scaled + 
## mm13:     wbf_cent_scaled:Treatment.order + wbf_cent_scaled:size_pc1_scaled + 
## mm13:     percLoad_cent_scaled:Treatment.order + Treatment.order:size_pc1_scaled
## mm6: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled + 
## mm6:     percLoad_cent_scaled + Treatment.order + size_pc1_scaled + 
## mm6:     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 | 
## mm6:     Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order + 
## mm6:     strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order + 
## mm6:     wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order + 
## mm6:     Treatment.order:size_pc1_scaled
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)   
## mm13 17 179.69 215.29 -72.844   145.69                            
## mm6  18 171.04 208.74 -67.518   135.04 10.652      1   0.001099 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mm14 <- update(mm6, .~. - strokeAmp2_scaled)
anova(mm6, mm14) # keep
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## mm14: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled + 
## mm14:     percLoad_cent_scaled + Treatment.order + size_pc1_scaled + 
## mm14:     wbf2_scaled + percLoad2_scaled + (1 | Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + 
## mm14:     strokeAmp_cent_scaled:Treatment.order + strokeAmp_cent_scaled:size_pc1_scaled + 
## mm14:     wbf_cent_scaled:Treatment.order + wbf_cent_scaled:size_pc1_scaled + 
## mm14:     percLoad_cent_scaled:Treatment.order + Treatment.order:size_pc1_scaled
## mm6: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled + 
## mm6:     percLoad_cent_scaled + Treatment.order + size_pc1_scaled + 
## mm6:     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 | 
## mm6:     Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order + 
## mm6:     strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order + 
## mm6:     wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order + 
## mm6:     Treatment.order:size_pc1_scaled
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)   
## mm14 17 179.11 214.71 -72.555   145.11                            
## mm6  18 171.04 208.74 -67.518   135.04 10.073      1   0.001504 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mm6)  # final model for paper
## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +  
##     percLoad_cent_scaled + Treatment.order + size_pc1_scaled +  
##     strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |  
##     Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order +  
##     strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order +  
##     wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order +  
##     Treatment.order:size_pc1_scaled
##    Data: bdta
## 
## REML criterion at convergence: 163.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6717 -0.3034 -0.0680  0.3528  2.2000 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 1.1046   1.0510  
##  Residual             0.2438   0.4937  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                            Estimate Std. Error t value
## (Intercept)                                10.33167    0.33224  31.097
## strokeAmp_cent_scaled                      -1.21537    0.36143  -3.363
## wbf_cent_scaled                             1.20177    0.26092   4.606
## percLoad_cent_scaled                        2.45005    0.38932   6.293
## Treatment.order2                            0.06881    0.15588   0.441
## size_pc1_scaled                             3.18673    0.29199  10.914
## strokeAmp2_scaled                           1.36399    0.49926   2.732
## wbf2_scaled                                -0.27572    0.24625  -1.120
## percLoad2_scaled                            1.12004    0.40566   2.761
## strokeAmp_cent_scaled:percLoad_cent_scaled -2.08682    0.58439  -3.571
## strokeAmp_cent_scaled:Treatment.order2      0.96049    0.44821   2.143
## strokeAmp_cent_scaled:size_pc1_scaled       0.19023    0.10033   1.896
## wbf_cent_scaled:Treatment.order2           -0.48672    0.26616  -1.829
## wbf_cent_scaled:size_pc1_scaled            -0.51630    0.21509  -2.400
## percLoad_cent_scaled:Treatment.order2      -0.98064    0.41651  -2.354
## Treatment.order2:size_pc1_scaled           -0.88736    0.23299  -3.809
## 
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it
# write output
summary(mm9)$coefficients  
##                                               Estimate Std. Error
## (Intercept)                                10.50939334  0.3719988
## strokeAmp_cent_scaled                      -0.78693814  0.4224585
## wbf_cent_scaled                             0.79953361  0.2846480
## percLoad_cent_scaled                        1.99184379  0.4507705
## Treatment.order2                            0.03758213  0.2003348
## size_pc1_scaled                             2.62560994  0.2567266
## strokeAmp2_scaled                           0.88098584  0.6199180
## wbf2_scaled                                -0.39713261  0.3031080
## percLoad2_scaled                            0.74256129  0.5006662
## strokeAmp_cent_scaled:percLoad_cent_scaled -1.30791510  0.7164999
## strokeAmp_cent_scaled:Treatment.order2      0.69012764  0.5682134
## strokeAmp_cent_scaled:size_pc1_scaled       0.17030338  0.1302434
## wbf_cent_scaled:Treatment.order2            0.18640327  0.2572330
## wbf_cent_scaled:size_pc1_scaled            -0.45047772  0.2380927
## percLoad_cent_scaled:Treatment.order2      -0.86679594  0.5231384
##                                               t value
## (Intercept)                                28.2511470
## strokeAmp_cent_scaled                      -1.8627583
## wbf_cent_scaled                             2.8088498
## percLoad_cent_scaled                        4.4187536
## Treatment.order2                            0.1875966
## size_pc1_scaled                            10.2272624
## strokeAmp2_scaled                           1.4211327
## wbf2_scaled                                -1.3102017
## percLoad2_scaled                            1.4831465
## strokeAmp_cent_scaled:percLoad_cent_scaled -1.8254227
## strokeAmp_cent_scaled:Treatment.order2      1.2145571
## strokeAmp_cent_scaled:size_pc1_scaled       1.3075783
## wbf_cent_scaled:Treatment.order2            0.7246475
## wbf_cent_scaled:size_pc1_scaled            -1.8920268
## percLoad_cent_scaled:Treatment.order2      -1.6569152
write.csv( summary(mm11)$coefficients, file = "resp_oth_Coefs_percLoad.csv" )

met rate with other predictors diagnostics

# qq plot
qqnorm(resid(mm6), main = "")
qqline(resid(mm6)) # good, but two outliers

# residual plot
plot(fitted(mm6), resid(mm6), xlab = "fitted", ylab = "residuals")
abline(0,0)

# QQPlot for group-level effects
qqnorm(ranef(mm6)$Bee.ID[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(mm6)$Bee.ID[[1]]) # looks good

# visualize model: 
sjp.lmer(mm6, type = 'fe')
## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.

met rate with other variable visualization

newdf <- data.frame(strokeAmp_cent_scaled = 0, 
                    wbf_cent_scaled = 0, 
                    percLoad_cent_scaled = seq(min(bdta$percLoad_cent_scaled), 
                                               max(bdta$percLoad_cent_scaled), 
                                               length.out = 100),
                    Treatment.order = factor(1), 
                    size_pc1_scaled = 0, 
                    strokeAmp2_scaled = 0, 
                    wbf2_scaled = 0
                    )
newdf$percLoad2_scaled <- scale(newdf$percLoad_cent_scaled^2)

preds1 <- predict(mm6, newdf, re.form = NA)

# plot of prediction for average sized bee (line)
# raw data plotted as points
ggplot(bdta, aes(x = percLoad_cent_scaled , y = av.resp..CO2.mL.hr.)) + 
     geom_point(alpha = 0.3) + 
     geom_line(data = newdf, aes(x = percLoad_cent_scaled, y = preds1)) + 
     labs(x = "Mass during flight (% of starved mass in g)")

# visualize metabolic rate vs. wbf
newdf <- data.frame(strokeAmp_cent_scaled = 0, 
                    wbf_cent_scaled = seq(min(bdta$wbf_cent_scaled), 
                                               max(bdta$wbf_cent_scaled), 
                                               length.out = 100), 
                    percLoad_cent_scaled = 0,
                    percLoad2_scaled = 0,
                    Treatment.order = factor(1), 
                    size_pc1_scaled = 0, 
                    strokeAmp2_scaled = 0
                    )
newdf$wbf2_scaled <- scale(newdf$wbf_cent_scaled^2)

preds1 <- predict(mm6, newdf, re.form = NA)

# plot of prediction for average sized bee (line)
# raw data plotted as points
ggplot(bdta, aes(x = wbf_cent_scaled , y = av.resp..CO2.mL.hr.)) + 
     geom_point(alpha = 0.3) + 
     geom_line(data = newdf, aes(x = wbf_cent_scaled, y = preds1)) + 
     labs(x = "Wing beat frequency (centered and scaled)")

# visualize metabolic rate vs. stroke
newdf <- data.frame(
                    wbf_cent_scaled = 0,
                    wbf2_scaled = 0,
                    percLoad_cent_scaled = 0,
                    percLoad2_scaled = 0,
                    Treatment.order = factor(1), 
                    size_pc1_scaled = 0, 
                    strokeAmp_cent_scaled = seq(min(bdta$strokeAmp_cent_scaled), 
                                               max(bdta$strokeAmp_cent_scaled), 
                                               length.out = 100)
                    )
newdf$strokeAmp2_scaled <- scale(newdf$strokeAmp_cent_scaled^2)

preds1 <- predict(mm6, newdf, re.form = NA)

# plot of prediction for average sized bee (line)
# raw data plotted as points
ggplot(bdta, aes(x = strokeAmp_cent_scaled , y = av.resp..CO2.mL.hr.)) + 
     geom_point(alpha = 0.3) + 
     geom_line(data = newdf, aes(x = strokeAmp_cent_scaled, y = preds1)) + 
     labs(x = "Stroke Amplitude (centered and scaled)")

This model is fundamentally different from the previous one. It’s saying, when we account for differences in wingbeat freq, stroke amplitude, and bee size, we sill see an effect of percent load on metabolic rate.

OR, we could say that holding all other variables constant, an increase in wingbeat frequency is associated with an increase in metabolic rate.